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' Andreev billiards are finite, arbitrarily-shaped, normal-state regions, surrounded by supercon- 

ductor. At energies below the superconducting gap, single-quasiparticle excitations are confined 
to the normal region and its vicinity, the essential mechanism for this confinement being Andreev 
reflection. This Paper develops and implements a theoretical framework for the investigation of the 
short-wave quantal properties of these single-quasiparticle excitations. The focus is primarily on 
\^ ' the relationship between the quasiparticle energy eigenvalue spectrum and the geometrical shape of 

the normal-state region, i.e., the question of spectral geometry in the novel setting of excitations 
confined by a superconducting pair-potential. Among the central results of this investigation are 
two semiclassical trace formulas for the density of states. The first, a lower-resolution formula, 
corresponds to the well-known quasiclassical approximation, conventionally invoked in settings in- 
volving superconductivity. The second, a higher-resolution formula, allows the density of states to 
be expressed in terms of: (i) An explicit formula for the level density, valid in the short-wave limit, 
for billiards of arbitrary shape and dimensionality. This level density depends on the billiard shape 
only through the set of stationary- length chords of the billiard and the curvature of the boundary at 
the endpoints of these chords; and (ii) Higher-resolution corrections to the level density, expressed 
as a sum over periodic orbits that creep around the billiard boundary. Owing to the fact that these 
creeping orbits are much longer than the stationary chords, one can, inter alia, "hear" the stationary 
chords of Andreev billiards. 
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I. INTRODUCTION AND OVERVIEW 



The purpose of this paper is to consider the quantal dynamics of elementary electron and hole quasiparticle ex- 
citations existing within and in the vicinity of a normal-state region of matter that is completely surrounded by an 

■ essentially infinite region of conventional superconductor. The entire system — normal-state region and superconduct- 
ing surround — may be envisaged as three-dimensional, although the approach that we shall be developing is applicable 

t-H , in any number of dimensions. Owing to the inability of the surrounding superconductor to support propagating quasi- 
particle excitations at sufficiently low energies, electron and hole quasiparticle excitations at such energies are bound 
to the normal-state region and its vicinity, and it is on the properties of such bound states that we shall be focusing 

■ our attention. 

The process responsible for the confinement of these excitations to the normal-state region and its vicinity is Andreev 
reflection Q from the surrounding superconducting condensate; we shall therefore refer to such structures, near to 
which quasiparticles are confined, as Andreev billiards. Andreev billiards were introduced and certain simple aspects 
I i of their classical dynamics were discussed in Ref. . A very brief account of the approach and results contained in 
the present Paper were reported in Ref. || . Certain quantum-mechanical properties of Andreev billiards were studied 
CI inRefs. f§. 

As we explore the quantal dynamics of quasiparticle excitations of Andreev billiards, our primary focus will be on 
the relationship between the quasiparticle energy eigenvalue spectrum and the geometrical shape of the normal-state 
^ ■ region, i.e., the question of spectral geometry in this novel setting of excitations confined by a superconducting pair- 
k>( \ potential. In the setting of conventional billiard systems || confinement is, by contrast, accomplished by an infinite 
(or occasionally finite) single-particle potential pqj . As mentioned above, we shall primarily be concerned with 
confined quasiparticle states, and therefore in energy eigenvalues lying within the quasiparticle gap of the surrounding 
superconductor (although our approach is also suited to the study of scattering states). The strategy that we shall 
develop is inspired by the beautiful work of Balian and Bloch, in which, inter alia, the eigenvalue spectrum of the 
Laplace operator was investigated for generically-shaped spatial regions and various types of boundary conditions ||-|| . 
The central theme of the work of Balian and Bloch is the relationship between the boundary shape, the type of 
boundary conditions, and the the eigenvalue spectrum. We shall refer to Refs. ||-[| respectively as BB-I, BB-II, and 
BB-III. 

As it is so central to the properties of Andreev billiards, let us pause to review the core qualitative features of 
the Andreev reflection process: to a high degree of accuracy it (i) interconverts electron and hole excitations; and 
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(ii) reverses the velocity of excitations. It is this latter, retro-reflective, character of the Andreev reflection process 
that endows Andreev billiards with dynamical characteristics quite distinct from those of conventional billiards, in 
which confinement is caused by specular reflection from a single-particle potential. 

We shall take as a model for the normal-state region of an Andreev billiard a Fermi gas, parametrized by the Fermi 
energy. Thus, in the normal-state region we shall be neglecting the effects of band structure, impurity scattering, and 
quasiparticle interactions. We shall account for the superconducting nature of the matter surrounding the normal- 
state region by asserting that there is a superconducting pair-potential A that varies discontinuously: inside the 
normal-state region A = 0; outside the normal-state region A takes on the constant value Ao 0). We shall 
refer to the surface on which A changes discontinuously as the shape of the Andreev billiard. Thus, we shall not 
be working self-consistently, but shall benefit from being in a position to develop an interface-scattering approach 
to the quasiparticle dynamics, in which we are able to focus on processes occurring at the interface. Hence, we can 
incorporate in a direct and natural manner the impact of the shape of the billiard (i.e. the shape of the interface) on 
the spectrum of energy eigenvalues of the confined electron-hole quasiparticles. 

We see four principal sources of motivation for the present work. First, Andreev billiards provide a novel setting 
for the exploration of spectrum-shape relationships, a branch of mathematics with a distinguished history |]l7| , |l8|| . 
The novelty is fed in by the Andreev reflection process occuring at the normal-to-superconductor boundary. Second, 
the usual spectral-geometric scenario [in which deviations of the density of modes from its large-system limit become 
appreciable as the wavelength become comparable to the characteristic linear size of the system] is not the whole 
story for the case of Andreev billiards. Instead, owing to the presence of a second, much larger, lengthscale, set 
by the difference between the momenta of incident electrons and the holes they become upon Andreev reflection 
(and vice versa). As this momenutm difference is small (on the scale of the Fermi momentum), the corresponding 
lengthscale is much larger that the Fermi wavelength. Through this new lengthscale the eigenvalue spectrum can be 
sensitive to the shape of the billiard even when the characteristic size of the billiard is much larger than the underlying 
wave length associated with quasiparticle motion. This relevance of the new lengthscale has long been appreciated, 
showing up, e.g., in the effectively one-dimensional settings of Tomasch jl9| and McMillan- Rowell [^0| oscillations 
in the tunneling density of states above the superconducting gap, and in de Gennes-Saint- James bound states |2l| l 
below the superconducting gap Third, as we shall see when we develop a trace formula for the (oscillatory part of 
the) density of quasiparticle eigenstates (DOS), there turns out to be a novel and useful separation in the scale of 
periods of the two dominating classes of (primitive, classical, periodic) trajectories that feature. As a consequence, the 
DOS will comprise: (i) a relatively smooth contribution due to retracings of geometrical chords across the billiard of 
stationary lengths, dressed by (ii) a more rapidly varying contribution arising from orbits located near the boundary 
and involving charge-preserving as well as charge-interconverting reflection processes. Thus, from the oscillatory part 
of the DOS one can "hear" aspects of the shape of the billiard such as the stationary values of the lengths of the 
chord ). We are not aware of any other spectral-geometric contexts that feature this type of information. Fourth, the 
quasiparticle energy eigenvalue spectrum, and its sensitivity to the shape of the billiard, should be experimentally 
accessible, e.g., via tunneling spectroscopy on hybrid superconducting/normal-state structures. The current state of 
microfabrication technology makes such experiments realizable |22f . 

We see the following as the principal results of the present work. First, we provide the machinery for computing the 
Green function for Andreev billiards of arbitrary shape in terms of a multiple scattering expansion that focuses on the 
influence of the billiard shape. Second, we implement this machinery to construct two semiclassical schemes (resulting 
in two semiclassical trace formulas) for computing the oscillatory component of the DOS. One, which we shall refer to 
as Scheme A, simply amounts to an elaboration (to billiards of arbitrary shape) of Andreev's approximation. Thus, it 
gives a DOS that takes the form of an integral over the chords of the normal-state region with an appropriate weight 
function. Hence, it realizes the intuitively natural notion that the chords, being the periodic orbits at the Andreev 
level, determine the energy eigenvalue spectrum, according to Bohr-Sommerfeld quantization conditions. The other, 
Scheme B, captures certain physical effects that are inaccessible to Scheme A, such as mesoscale oscillations in the 
DOS. 

This paper is organized as follows. Following the present introduction and overview we define, in Sec. [n|, Andreev 
billiards and present the corresponding Bogoliubov-De Gennes (BDG) eigenproblem. In this section we also intro- 
duce the Green function for the BDG eigenproblem and provide its connection to the DOS, and review the standard 
quasiclassical approach to the BDG eigenproblem, due to Andreev. In Sec. [II we formulate the computation of the 
Green function in terms of an expansion in which the basic processes are reflection from and transmission through the 
interface separating the normal(N) and superconducting(S) regions. We then make a physically-motivated reorgani- 
zation of this expansion, which will later allow us to integrate out states in the S region, thus obtaining a description 
solely in terms of states within the billiard. (Section [nj is the technical basis of the Paper; however, it may safely 
be omitted by readers wishing to focus on results rather than methods.) After this exact reformulation we proceed, 



in Sec. IV, to integrate out the states in the S region within an approximation scheme valid for short waves. In 
this way, we obtain an effective Multiple Reflection (rather than Scattering) Expansion. We continue, in Sec. |v|, by 
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following certain approximation strategies that allow us to compute the DOS at various levels of energ y- re solution, 
in each case obtaining the corresponding trace formula. This section is the heart of the Paper. In Sec. VI we make 
some concluding remarks and hint at some possible applications of the ideas we have presented. We have relegated 
to appendices some background material, including a derivation of the BDG wave equation and the rudiments of 
boundary integral methods, as well as some technical and parenthetical passages. 



II. ANDREEV BILLIARDS 



A. Idealization of the physical system 



The physical system of interest in the present Paper is an Andrccv billiard of arbitrary shape, i.e., a normal-state 
region embedded inside an infinite superconducting region, as depicted in Fig. ^. Following Gor'kov's mean-field 
approach to superconductivity [ p3"|]24f| , we describe the system by the (variable particle-number) Hamiltonian H, 
given by 



H =Hj ^^(x){-|^V 2 - A1 |v' Q (x)+ J d d x {A(x)> + (x)^_(x)+A(x)4(x)V>!(x)} 



(2.1) 



Here, m is the effective electron mass, fi is the electron chemical potential (which we take to be uniform throughout the 
system), A(x) is a given superconducting pair-potential which characterizes the superconducting condensate, V\L( X ) 
and ^ Q (x) are creation and annihilation field operators for electron quasiparticles having position x and spin-projection 
a = ±, and d is the dimension of space (typically two or three). We shall not go beyond the picture of quasiparticle 
excitations propagating in the presence of a superconducting condensate implied by this description. Within the 



superconducting 
region (S) 




FIG. 1. Two-dimensional example of an Andreev billiard, showing a normal region (N) surrounded by a superconducting 
region (R). In this example the billiard is convex. 

Gor'kov description of the consequences of the electron-electron interaction, any (Heisenberg-representation) excited 
state that is arrived at by the addition of a single spin-up electron quasiparticle to the (Heisenberg-representation) 
ground state |<I>o) evolves into a coherent superposition of such a state and a stat e arrived at by the removal of a 
spin-down electron quasiparticle from the ground state p5| ], the Hamiltonian ( |2.1[ ) maintaining the system in this 
sector of Fock space. Thus, it is adequate to address states of the form 

1*0 = f d d x (u(x)$(x)+w(x)Vi(x)) |3o>, (2.2) 
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which are described by the two-component complex-valued amplitudes (u(x), u(x)J , i.e., the family of one-quasiparticle 
excited states. 

To derive the equation of motion for these one-quasiparticle states (i.e. the so-called time-dependent BDG equa- 
tion P|| ), we follow Andreev jjj and suppose that the system is in some (Heisenberg- representation) one-quasiparticle 
excited state Then the amplitude u(pc,t) for finding the state at time t to be the ground state with an up-spin 

electron quasiparticle added at position x (i.e. the Heisenberg-representation state ^_J_(x, t)|$o)) is given by 

($ 1 V>+(x, <)|$!>. (2.3a) 

Similarly, the amplitude w(x, t) for finding the state at time t to be the ground state with a down spin electron 
quasiparticle removed at position x (i.e. the Heisenberg-representation state ^-(x, t)|$ )) is given by 

($o|V>l(x,i)|$i>. (2.3b) 

Here, ip + (x,t) = e lHt / h ?/> + (x) c~ lHt / h and "0-( x , t) = Q tHt / h -0^ ( x ) e ~ tHt / h ar6j respectively, Heisenberg- 
representation field operators. Thus, the wave functions u(x, t) and v(x, t) serve as amplitudes for the present 
up-spin electron quasiparticle and the absent down-spin electron (i.e. up-spin hole) quasiparticle. Then, by virtue 
of the He isenberg equation of motion for the field operators (see, e.g., Ref. pj], Sec. 6), together with the Hamilto- 



nian (2.1), it is a straightforward exercise in computing commutators of field operators to show that the amplitudes 
w(x, t) and w(x, t) evolve according to the appropriate time-dependent Schrodinger equation, i.e., the time-dependent 
BDG equation: 

a /-I^v 2 - M A( X ) 




Analysis of this equation via the separation of the time variable, appropriate when there is no external time- 
dependence, leads to the BDG eigenproblem 





= E n \ , (2.5) 



A*(x) 

where n is a (collective) index for all quantum numbers and {E n } and { (it n (x), i; ra (x)) } are the corresponding energy 
eigenvalues and (two-component) eigenfunctions. 

In general, A may vary spatially. However, we shall consider the situation in which deep inside the superconductor A 
goes to a constant value Ao, whereas throughout the normal metal it vanishes. In the intermediate region (i.e. within 
a superconducting coherence length outside of the billiard boundary) A is suppressed to a value smaller than Ao, 
and falls to zero as the N region is entered. We shall ignore the effects of this suppression and assume that A varies 
discontinuously between to Ao at a surface, which we refer to as the billiard boundary and denote by dV, that 
divides the system into two homogeneous regions, the billard interior (denoted V) and the billiard exterior (denoted 
V). For the sake of simplicity, we further assume that there are no metallurgical differences between the normal-state 
and superconducting regions, inasmuch as the only difference between them is the value of the pair potential (the 
effective mass, e.g., being common). 

To ease the notation we shall adopt units in which ft 2 /2m = 1. To recover results in terms of the original physical 
units, one multiplies the three variables having the dimensions of energy [viz. /i, A(r) and E] by the factor (2m/S 2 ). 

The BDG eigenproblem plays the same role for Andreev billiards that the Schrodinger eigenproblem plays for 
conventional billiards. If a conventional billiard is surrounded by a region in which the single-particle potential is 
infinite then we refer to the billiard as a hard billiard, and the Schrodinger equation outside the billiard is replaced 
by the homogeneous Dirichlct (i.e. vanishing) boundary condition on the Schrodinger eigenfunction, this boundary 
condition leading to the quantization of the eigenvalue spectrum. If, on the other hand, a conventional billiard is 
surrounded by a region in which the single-particle potential is finite then we refer to the billiard as a soft billiard, 
and the solution to the Schrodinger equation outside must be matched on to the the solution of the Schrodinger 
equation inside, this matching leading to the quantization of the eigenvalue spectrum. The Andreev billiard is, 
therefore, analogous to a soft Schrodinger billiard; its hard limit seems difficult to realize because, at least in known 
superconductors, the pair potential is far smaller than the chemical potential. 



The eigenproblem for Andreev billiards, then, is given by Eq. (2.5) with A(x) = for x inside the billiard, and 



A(x) = A for x outside the billiard. Thus, we are faced with the task of addressing the BDG eigenproblem for the 
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case in which the system comprises two spatially homogeneous regions (one normal, one superconducting) that meet 
at a closed surface. Owing to the spatial homogeneity of these regions, the general solution of the BDG equation can 
readily be obtained in each region. The quantization of the eigenvalue spectrum results from the matching of the 
solutions and their normal derivatives across this surface, together with the confinement of the eigenfunctions to the 
vicinity of the Andeev billiard. The resulting spectrum depends, therefore, on the shape of this surface. Exploring 
this dependence is the central aim of the present work. It would be straightforward to extend the present framework 
to handle issues such as Josephson coupling between superconducting regions, scattering from Andreev billiards, etc. 



B. Green function and density of states 

The spectrum of eigenvalues of the BDG wave equation {E n } is assembled into the DOS p(E), which is defined as 

p(E) = Y,S(E-E n ). (2.6) 

n 

As is commonly the case, it is convenient to approach p(E) via a Green function for the BDG wave equation, which 
we now introduce. 



1. Green function for the BDG equation 

The (2x2 matrix) Green function G(x, x'; z) for the BDG wave equation is defined by the following matrix partial 
differential equation: 



V — n — z A(x) 
A*(x) V 2 +p-z 



G(x,x';z) = I5(x-x'), (2.7) 



where z is the complex energy and I is the (2x 2) identity matrix. Under the far-field boundary condition G(x, x'; z) — > 
as | x — x' | — > oo the Green function G(x, x' ; z) is unique |^7j . For the case of the Andreev billiard we have A (x) = 
for x in V and A(x) = Ao for x in V. 

It is useful to express G(x, x'; z) in terms of the eigenfunctions (u„, v n ) of the BDG Hamiltonian, i.e., 



w„(x)u*(x') u„(x)«*(x') 

«n(x)<(x') V n (x)v*(x') 



(2.8) 



In order to see that this form does indeed satisfy Eq. (2.7), one may substitute this expression into Eq. (2.7) and 
make use of Eq. ( [2~5| ) and the completeness of the eigenfunctions, i.e., 

/u„(x)<(x') u„(x)u*(x')\ 

„ Un(x)<(x') t,„(x)<(x'W 



2. Connection between BDG Green function and density of states 



We now follow the standard practice of expressing a DOS in terms of the corresponding Green function by making 
use of the identity 



— Im lim 

7T E„ — E — it 



= S(E n -E). 



(2.10) 



Together with Eqs. (2J3) and (2J:), this allows us to see that 

p{E) = \ fd d x {ImTrG(x,x';£ + !£ )} x , =x = J d*x £ Ilm^^ ~ H ' '"' X) 



E-ie 
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Here, TrG d enote s the (matrix) trace over the diagonal components of G. 



Equation ( 2.11 ) is an expression for p(E) in terms of G(x, x.';E + ie), and any approximation to G(x, x.';E + ie) 



thus furnishes an approximation to p(E). However, being a sum of delta functions, p{E) is not a smooth function 
and can, therefore, be extremely awkward to approximate. In order to find approximations to p{E) it is preferable 
to seek a continuous function that carries essentially the same information as it. One candidate is the smoothed DOS 
p 1 (E), defined via 



Pl (E) 



dE' f(E-E'; 7 )p(E'), 



(2.11) 



where f(E — E';j) is some smoothing function and 7 is the (real) smoothing width. We remark, parenthetically, 
that any DOS derived from experiment will be smoothed to some extent, depending on the resolving power of the 
apparatus and/or the lifetime of the single-particle excitations. 

There are several possible choices for the smoothing fuction f(E — E';j), including, e.g., Lorentzian, Gaussian and 
logarithmic-Gaussian. It is also possible to define a continuous integral transform of p(E) that is, itself, a physical 
property, and for which we can derive some approximation. For example, under the Lambert transform of p{E), which 
exchanges energy for temperature, p(E) is transformed into the total (equilibrium internal) energy as a function of 
temperature j28| . Which smoothing procedure is the best choice depends on the method used to approximate p(E). 
For Green-function-based methods, such as the one we shall adopt, the Lorentzian smoothing function, 



f(E-E'^) = - 



7 



TT (E - E'f + 7 2 ' 



(2.12) 



is the most appropriate, for reasons that should become clear below. 

We now give the analogue of Eq. ( 2.11 ) for relating the Lorentzian-smoothed D OS p -y{E) and the Green function at 
complex energy G(x, x'; z). Following essentially the procedure that lead to Eq. (2.11), we have the following identity 
for arbitrary (real) 7: 



pJE) = - [ d d xlmTiG(x,x';E + ij) = [ d d x V - Im 

7T J x'=x J *—? 7T 



|tt n (x)| 2 + |f w (x)|^ 
E n - E - i 1 



E 1 



TT (E - E n f + 7 2 



dE 



, 1 7 
it (E- E') 2 +7 



"£S(E'-E n ). 



(2.13) 



Naturally, in the limit 7 — > + , p-y(E) passes to p(E). 



3. Fundamental Green function for a homogeneous normal region 

In this section we derive the fundamental (i.e. homogeneous) normal-state Green function, the first of two Green 
functions that are central to the construction of the Green function for an Andreev billiard. Along the way, we 
introduce some convenient notation that we shall subsequently make use of. This fundamental Green function is the 
translationally-invariant solution G N (x, x';z) of the equation 

(H — zT) G (x, x'; z ) = 5^(x-x')I, (2.14a) 

H=(*V -p 2 + ,)" ( ^ 2 + ^ ),T3 ' (2 ' 14b) 

where <t 12 ,3 are the three Pauli matrices, together with the boundary condition that the Green function vanishes at 
infinity: 

G N (x,x';z)^0 as |x-x'|->oo. (2.15) 
The solution for this Green function is given by the (spatial) matrix element of the operator (zl — H) _1 , i.e., 

G N (x,x';,) = <x|^|x') = W ^^W 

= (zl - (Vl + p) <r 3 ) -L ( (x |(z p 2 + M )-V> + (x|(* + P 2 /i)-V>) , (2.16) 



G 



which is, of course, a matrix in electron/hole space. The two terms inside the parantheses are quite familiar: they 
are the usual Green functions of the Helmholtz wave equation, and can be evaluated in the standard way. In three 
dimensions, e.g., one has 

i p±ife±|x— x 

(x|(z T P 2 ± M)-V> = ±^ |x „ x1 , (2-17) 

The symbols k±, which will be used throughout this Paper, denote particle and hole wave numbers, and are given by 
the expressions 



k± = VfT±z , (2.18) 

where it is understood that the roots having positive real parts are the ones that are adopted. Then, for the three- 
dimensional case, one arrives at the the fundamental Green function for the normal state: 

N , zI-(V 2 + M )T3 1 e-|M*-^l \ _ 1 ( e *+l«-*'l \ , 9 1Q v 

lX ' X ' Zj ~ 2z 47r|jx-x'| |X-X'| J ~ 4tt|x-x'| ^ _ e -ifc-|x-x'| ) ■ V'M) 

This Green function has the expected form: an outgoing spherical wave in the particle component (i.e. the particle 
Green function) and an incoming spherical wave in the hole component (i.e. the hole Green function), the latter 
representing an outgoing hole. Moreover, as we expect in the normal state, there is no mixing between the particle 
and hole components, the off-diagonal elements being zero. Next, we introduce some convenient notation for the 
components of this fundamental Green function, viz., 

p±ife± |x— x'| 

ff± ( x , x ') = - j-, (2.20) 



in terms of which G N becomes 



4- Fundamental Green function for a homogeneous superconducting region 

The fundamental Green function for a homogeneous superconducuting region is the translationally invariant solution 
G s (x, x'; z) of the equation: 

{zl + (V 2 + n) (T 3 + Ann + A 2 <r 2 } G s (x - x') = I«5 (d) (x - x'), (2.22) 

together with the boundary condition that it vanishes as |x — x'| goes to infinity. Here Ai and A 2 represent the 
(constant) real and imaginary parts of the complex pair-potential. Note that we shall henceforth take the complex 
energy z to be the real energy E. We shall be concerned with situations in which the quasiparticles are bound to the 
billiard and shall, threfore, assume that \E\ < A, where the magnitude A of the pair potential obeys A 2 = A 2 + A 2 ,. 
In the present work, it is convenient to choose a specific gauge, which we do by setting Ai = and A 2 = A. However, 
for extensions of the present work to settings, such as SNS junctions, in which there are physical implications of phase 
differences it is necessary to ksow the Green function for arbitrary gauges, and it is therefore for this case that we 
provide the Green function. 



To obtain the Green function one first formally inverts Eq. fl2.22|) to obtain 



G s (x-x') = {£I+(V 2 +//)o-3 + Ai<7i + A 2 o- 2 } ^(x-x'). (2.23) 

By manipulating this equation so as to separate the matrix and partial-differential aspects of the inversion operation 
one then arrives at 

G s (x-x') = {^I+(V 2 +^)o-3 + A 1 cr 1 + A 2C r 2 } _1 5 3 (x-x') (2.24a) 
= {EI - (V 2 + /i)<7 3 - A lCTl - A 2 er 2 } {E 2 - A 2 - (V 2 + fi) 2 }^ <5 (d) (x - x') (2.24b) 
= {EI- (V 2 +h)(t 3 - Aio-i - A 2 cr 2 } 



2WA 2 - E 2 



i\/A 2 -£ 2 ) 1 - (v 2 + fi - iy/A 2 - £ 2 ) J j £(«9( x _ x '). (2.24c) 
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We identify the wave vectors fc+ and k^_ for the electron-like and hole-like components, respectively, as follows: 



k% = \J (i±i\/A 2 - E 2 , Refc|>0. 



(2.25) 



As the partial differential operators (v 2 + (4) 2 ) 

are partial differential operators of the Helmholtz wave equati< 
for wave vectors of length k±, their inversion is well known, 



ffl(x-x') 



V 2 



k±) ) ^ 3 (x-x') 



g±ifc|_ |x— x'| 

47r|x - x'l 



(2.26) 



where the final form on the right hand side holds for the three-dimensi onal case, and the minus sign in the exponent 
for holes ensures the proper decay at large distances. Upon using Eq. (2.26), G s becomes 



G s (x-x') = {EI- (V 2 +n)<73- Ai<n - A 2 <r 2 } 



1 



E 



2 [WA 2 -E 2 WA 2 - E 2 



2WA 2 - E 2 
A 2 



(4(x-x')-<£(x-x')) 



iVA 2 - E 2 



<r 2 |(^(x-x')-.9 S (x-x')) 



+ -<X 3 (4(X-X0+£(X-X')) 



(2.27a) 



(2.27b) 



To get to the second line, the action of (V 2 + /i) on g± is calculated with the help of Eq. (2.26). 



C. Andreev's quasiclassical aproximation scheme 

We now review the conventional approximation scheme for studying N-S hybrid systems, such as Andreev billiards. 
This approximation scheme was put forth by Andreev jjj, and we shall refer to it as Andreev's approximation. It 
takes advantage of the fact that for typical N-S systems the dimensionless parameter Af[i is much less than unity. 
The scheme consists of the separation of rapid and slow oscillations in the wavefunction. It can be motivated in the 
following way: consider an arbitrary N-S structure (i.e. arbitrary A(x) <C [i). Suppose that the electronic properties 
of the system are probed during a short time At, so that there is insufficient energy resolution around E = to resolve 
A(x). Then no measurement obtained via this probe is capable of distinguishing between the original system and a 
system described by the same BDG equation but with A(x) set to zero. An energy resolution of AE around E = 
corresponds to a momentum resolution of Ap — m AE/kp which, via the Heisenberg uncertainty relation, corresponds 
to a spatial resolution of Ax — h 2 k-p/m AE. Therefore, in order to resolve any effects of the pair-potential on the 
electronic states, the system must be probed on lengthscales larger than £ = fi^p/mA, i.e., the superconducting 
coherence length. Thus, the single-particle wavefunctions must have the form of plane- waves at the Fermi momentum 
with an envelope that varies on the lengthscale £, i.e., 



"(x) \ = e *fc F „.x ( u(x) 
u(xW U(x) 



(2.28) 



where the unit vector n defines the orientation of the wavevector of the plane wave and u and v are the slowly varying 
envelope amplitudes. By substituting this form into the BDG equation (2.5) one obtains 



- (V 2 + 2ik F n ■ V) 
A*(x) 



A(x) 
(V 2 + 2ik F n ■ V) 



As u and v vary on the lengthscale £, one has that 

|V 2 ?2| 



|fcFn • Vu| 



O 



w(x) 
v(x) 



< 1 



E 



u(x) 

e(x) 



(2.29) 



(2.30) 



Thus, to leading order in A/// it is permissible to ignore the V 2 term in Eq. ( 2.29 ), a singular approximation because it 
involves changing the order of the system of partial differential equations. This approximation to the BDG Hamiltonian 
is Andreev's approximation, and leads to the Andreev eigenproblem 
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-2ik F n-V A(x) 
A*(x) 2ik F n-V 



u(x) 



(2.31) 



Let us now apply Eq. ( 2.31 ) to Andreev billiards. First, it is useful to express the variable x in terms of b (i.e. the 
impact parameter or, equivalently, the transverse parameter), and s (i.e. the longitudinal parameter) such that 

x = b + ns. (2.32) 

Notice that b represents the trans verse degree(s) of freedom whereas s represents the longitudinal degree of freedom 
of the excitation. As, in Eq. ( 2.31 ), there is no differential operator acting on the variable b, the wavefunctions take 
the form 



u(h, s) 
v(h, s) 



J(b-b ) 



u(s;b ) 
v(s;b Q ) 



(2.33) 



where bo can be interpreted as the transverse quantum number. By substituting this form into Eq. (2.31) one obtains 

(2.34) 



-2ik F d/ds A(b ,s) 
A*(b ,s) 2ik F d/ds 



u{s\h ) \ _ E 



v(s;b ) 



u(s;b ) 
v(s;h ) 



Thus, one has reduced the partial differential eigenvalue equation (2.5) to a family of approximate ordinary differential 
eigenvalue equations parametrized by (n, bo). 

Now, for Andreev billiards A is real and piecewise constant. In this case, the solution of Eq. ( 2.31| ) proceeds as 
follows. The parameters (n, b ) fix a line, which determines A(b , s), in the sense that A(b , s) == A(b + ns). As 



A(bo, s) is piecewise constant, the task of solving Eq. (2.34) is straightforward. The quantization condition depends 

le line (specified by (n, b )) lying inside N]: 



on the chord length ^(n, bo) [i.e. the length of the part of t 

E£(n,b )/k F - 2ip = 2?rm 



(2.35) 



where m = 0, ±1, ±2, . . . and ip = cos _1 (_E/A). Thus, for each chord there is a ladder of energy eigenvalues. In order 
to obtain the DOS one must first obtain the DOS for a single chord, and then sum over all chords. However, it is 
not a priori completely clear w hat w eight should be assigned to each chord in performing the continuous summation. 
In fact, as we shall see in Sec. VA, at the Andreev level of approximation the DOS for a (convex, d-dimensional) 
Andreev billiard is given by 



p(E) 



°° hd-2 f / F 

£ 2(2^ ydndb%,b )^-£(n I bo)-2^-2 7 



(2.36) 



the two (d — l)-dimensional integrations, one over n and one over b, implement the summation over chords. 



III. MULTIPLE SCATTERING EXPANSION 



We now construct a M ultiple Scattering Expansion capable of yielding the Green function G(x, x') associated with 
the BDG equation (2/7), in settings in which A(x) is piecewise constant (i.e. takes on certain constant values in 
various regions of space that are delineated by d— 1-dimensional surfaces). Although the construction is applicable to 
a wider range of settings, we shall have in mind the application to an Andreev billiard, for which A(x) vanishes inside 
the billiard and has a constant nonzero value Ao outside it. The spirit of our approach very much parallels that of 
BB-I, although there there are significant differences arising from (i) the form of the eigenproblem (BDG rather than 
Laplace), and (ii) our need to employ matching (rather than boundary) conditions. 

In this section we set up the general formalism for the exact BDG Green function for the case in which the system 
is divided, via a change in the pair potential, into at least two distinct regions. As we shall see, the exact Green 
function can be expressed as a sum, generated by an interation scheme, over all possible scatterings from the boundary 
dividing these regions. 

Readers not familiar with the elements of potential theory that we shall be using (which are sometimes referred to 
as boundary integral techniques) may wish to pause to read App. in which we give a self-contained introduction to 
this subject and develop the elaborations necessary for application to the BDG eigenproblem. Specifically, we shall 
need to handle two-component wave functions in the context of matching (rather than boundary) conditions. 
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A. Matching conditions and boundary integral equations for the BDG Green function 



We now introduce a convenient parametrization of the Green function G(x, x'). We do this by decomposing G(x, x') 
into a particular integral and a complementary function. The particular integral, which yields the delta function under 
the action of the BDG operator, is built from the fundamental Green functions G N (x — x') and G (x — x'). The 
complementary function, which obeys the BDG wave equation [i.e. the homogeneous version of Eq. ( |2.7| , is specified 
in terms of as-yet undetermined single and double layers /i. n (a, x'), f OI (a,x'), /x IO (c*,x'), and i/°°(a,x'). Thus we 
write 



G(x,x') = < 



' G u (x,x') = G N (x-x') + J dv da a d a G N (x- a) /x u (a,x'), if x e V and x' e V; 

G oi (x,x') = f gv da a G s (x-a)u oi (a,x'), if x £ V and x' £ V; 

G io (x,x') = J 9v da a d a G N (x- a) /x io (a, x'), if x e V and x' e V; 

. G°°(x, x') = G s (x - x') + J gv da a G s (x - a) v°°{a, x'), if x £ V and x' e V. 



(3.1) 



We are employing the notation d a f(ct) to indicate the value of the component of the gradient of /(x) with respect 
to x, evaluated at the surface point a, directed along the inward normal direction at a. Similarly, da a indicates the 
(scalar) surface element at the point a. We shall find it useful to decorate G(x,x') with labels (such as, e.g., ii and 
io) according to the region of space in which the variables x and x' lie. For example, if x 6 V and x' <E V then we 
write G IO (x, x') for G(x, x'), conveying the idea that x lies inside V whereas x' lies outside V. 

Next, we fo cus on the surface dV across which the pair-potential is discontinuous. From the Green function 
equation ( |2.7| ) it is evident that both G(x, x') and its normal derivative n • V x G(x, x') are continuous as x varies 
across dV. What this means is that 



lim G(x,x')= Jim G(x,x'); 
lim ng V x G(x, x') = _lim na Va;G(x, x'). 



(3.2a) 
(3.2b) 



We apply this pair of matching conditions across dV to the parametrizations (3d), once fo r x' € V an d once for 
x' £ V. To evaluate the nec essary limits we make use of the using the continuity conditions (A31a) and ( A30b| ), as 
well as the jump conditions ( A30a ) and ( A31b ), thus arriving at the following conditions on the single layers u m and 
u°° and double layers fi u and ^t 10 : 





-G N (/3- 


-x')- 


/ da Q 

JdV 


9 Q G N (/3 


- a) fj. u (a 


x') + 


/ d<r a 

JdV 


G s (/3- 


a)^ oi (a,x'); 


(3.3a) 


i CT3 ^ oi (7,x') = 


9 7 G N ( 7 - 


x') + 


/ do a 

JdV 


a+a Q G N ( 7 - 


- a) fi n (a 


x')- 


/ da a 

JdV 


a 7 G s ( 7 - 


a)i/ oi (a,x'); 


(3.3b) 


i CT 3M io (/3,xO = 


G s ( 7 - 


x') + 


\ da a 

JdV 


G s (/3- 


a)v 00 {a, 


x')- 


I da a 

JdV 


d a G N ((3- 


a)/i io (a,x'); 


(3.3c) 


i < T 3 ^ 00 (7,x')- 


-9 7 G S ( 7 - 


x')- 


I da a 

JdV 


9 7 G S ( 7 - 




x') + 


/ da a 

JdV 


d+d a G N ( 7 - 


^a)/x io (a,x') 


(3.3d) 



These four matching conditions constitute a system of coupled integral equations for the single layers v m and v°° 
and the double layers fi n and fi la . It is convenient to collect together the single and double layers into a 4 x 4 matrix 
M(7,x'), defined via 



M(7,x') 

In terms of M(7,x') the system becomes 



'/x ii (7,x') /x io ( 7 ,x')' 
1/^(7, x') ^°°(7,x') 



M( 7 , x') = 2M°(7 - x') + 2 / da a G(-y, a) M(a, x'), 

JdV 



(3.4) 



(3.5a) 



where the inhomogeneity M° and the kernel G are defined via 
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l (7-x') = 



-<t 3 G n ( 7 -x') <x 3 G s ( 7 -x') 
<t 3 g> 7 G n ( 7 -x') -<r 3 d 7 G s ( 7 -x') 



(3.5b) 



G(7, a) 



-<r 3 d a G N ( 7 -a) <t 3 G s ( 7 -ck) 
cr 3 a 7 9 Q G N (7-a) -cr 3 9 7 G s (7-a) 



(3.5c) 



Equation (3.5a) is one of the central elements of this Paper. Its output (/x 11 etc.) can be fed into Eq. (3.1) to obtain 
the ingredients (G n , etc.) of G In terms of our 4x4 notation, this connection becomes 



'G u (x,x') G io (x,x')' 
G oi (x,x') G°°(x,x') 



'G N (x-x') 




G s (x-x') 



9 7 G N (x- 7 ) 



G s (x- 7 ) 



M(7,x') 



(3.6) 



When the iterative solution of Eq. (3.5a) is fed into Eq. (3.6) for the Green function, the resulting expansion is 




FIG. 2. Typical term in the MSE for the Green function. Lines running internally (externally) to the billiard represent 
homogeneous-region normal (superconducting) Green functions G N (G s ); each point on the boundary (a, /3, etc.) at which a 
scattering event occurs is to be integrated over the complete boundary. 

called a Multiple Scattering Expansion (MSE) for the Green function. A typical term in this expansion is shown 
diagrammatically in Fig. The physical content of the MSE is this: the iterations generate terms that correct 
the free Green function by accounting for multiple scatterings from the superconducting condensate surrounding the 
billiard. However, the expansion is not simply a perturbation expansion in powers of the pair-potential; instead, 
terms involving n transmissions (i.e. terms with G N followed by G s and vice versa) account nonperturbatively for 
all Feynman trajectories that traverse the boundary n times, thus spending intervals in the superconducting region. 
(The reason such a re-organization of the simple perturbation expansion in powers of the pair-potential is possible is 
that one knows fully the fundamental Green functions that describe propagation in homogeneous N or S regions.) 



B. Reorganized multiple scattering expansion 

The multiple scat tering expansion for the Green function was constructed by iterating a four-by-four matrix integral 
equation, viz., Eq. ( [3.5a ). In fact, this 4x4 structure consists of two substructures: (i) the electron-hole structure, 



which is an essential ingredient in Andreev billiards. (In fact it is an essential ingredient for any system involving 
superconductivity.) (ii) the inside-outside structure: this structure is specific to Andreev billiards, as they consist of 
two distinct regions separated by a boundary. As for the inside-outside structure, it is possible to diagonalize this (as 
we shall soon show), and thus to obtain a 2 x 2 matrix integral equation, whose iteration produces exactly the same 
MSE as the 4x4 matrix integral equation does. In this way, we reduce the problem of determining the full-space 
Green function to an effective, but nonetheless exact interior (or, if one wishes, exterior) problem. (Indeed, in the 
following section we shall obtain effective boundary conditions for the interior Green function problem, the solution 
of which coincides with the that for the full-space problem.) 

The advantages of this formulation are two-fold: First it allows us to easily obtain the (Feynman) rules for evaluating 
a generic term in the MSE. Second, as we shall see when we develop approximation schemes for the Green function, 
it is especially well suited for the task of integrating out the outside-propagation processes and, thus, obtaining an 
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effective Multiple Reflection Expansion. In this way, the physics of Andreev reflection as well as the corrections 
associated with charge-preserving reflection will become evident. 

To make this reorganization, let us introduce the operators ID) (diagonal) and O (off-diagonal), defined as follows by 
their action on 4-component functions: 



dap 



<0*(7) = 2 dap 



-<T 3 5 Q G N ( 7 -/3) 

_ CT3 a 7 G s (7-/3) 

< t 3 G s (7-/3) > 



ff3 5 7 5 Q G N (7-/3) 







*(/3). 



(3.7a) 



(3.7b) 



These operators constitute the (two-by-two block) diagonal and off- diagonal operator elements of the kernel G in 
the integral equation ( p. 5a ), the iterati ve sol ution of which generates the MSE. Note that we write the four-by-four 
identity as I. In terms of ID and O Eq. (3.5a) can be written symbolically (i.e. using an obvious condensed notation) 
as 



M = 2M° + 2 (I 



To obtain the reorganized MSE, first we rewrite Eq. (3.£) as 

(I - 2©) M = 2M° - 
and then, by inverting the operator (I — 2D), we obtain 

M = (I - 2D) -1 2M° + (I — 2D)" 1 



(3.8) 
(3.9) 
(3.10) 



in which the kernel has become (block) off-diagonal. Next, in order to obtain a (block) diagonal structure, we iterate 
this equation once, thus arriving at 



M = (I - 2D) 1 2M° + (I - 2D) 1 20 (I — 2D) 1 2M° + (I - 2D) 1 20 (I — 2D) 1 
At this stage it is useful to define the kernel K obeying 

(I - 2D) IK = I . 



This kernel is block diagonal, 



and the diagonal (two-by-two) blocks obey 



K n 

K°° 



K 11 + 2er 3 dG K il = I, 



K° 



2<r 3 <5G S K°° = I. 



(3.11) 
(3.12) 

(3.13) 



(3.14a) 
(3.14b) 



Here and elsewhere we shall use S (resp. d) without a subscript to indicate an inward normal derivative with respect 
to the first (resp. second) argument of the succeeding Green function, i.e., 



8G(a,P)=d a G(a,P), 
dG(a,/3) = df}G{a,p). 



In terms of K, Eq. (3.11), when pre- multiplied by (I — 2D), becomes 

(I - 2D) M = 2M° + 4G> IK M° + 40 IK O M, 
the upper-left (two-by-two) block of which can be rearranged to read 

H" = -2er 3 G N + 4<x 3 G s K°° <r 3 <5G N + {-2cr 3 <9G N + 4cr 3 G s K°° cr 3 86G N } jU U . 



(3.15a) 
(3.15b) 



(3.16) 



(3.17) 



What we have accomplished via these transformations is the construction of a closed (2x2) system of integral equations 
for the boundary layer fi u , which is all that is needed to complete the computation of the Green function G n . The 
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virtue of this transformation is that it facilitates the subsequent elimination (via the integrating out of processes 
involving virtual propagation in the superconducting region) of the states located in the superconducting region. This 
elimination can now be made straightfor wardl y, owing to the fact that all superconducting Green functions are now 



conveniently located in the kernel in Eq. ( 3.17 



The reorganization just described also allows one to identify the following rules for the construction of all possible 
contributions at any order n (= 1, 2,3,.. .) to the inside-to-inside Green function G 11 : 

1. Write down all possible permutations of G N and G s (having a total of n + 1 Green functions), subject to the 
restriction that the first and last Green functions are G N . 

2. Associate to each permutation a factor (— 1) JN+1 , where in is the number of G N factors. 

3. Furnish all G N factors (except the last) with normal derivatives acting on their second arguments; all G s factors 
carry no normal derivatives. 

4. Furnish any Green function factor that follows a G factor with an additional normal derivative acting on its 
first argument. 

5. Insert a Pauli-matrix factor <r 3 before every G N and G except the first. 

In this way one can construct G". An example is provided by the process depicted in Fig. || for which the corresponding 
amplitude is 

f do~ a dap da 1 da s <9G N (x, a) <x 3 dG N (a, (3) <x 3 <9G N (/3, 7) <r 3 G s ( 7 , S) <x 3 <SG N ((5, x'). (3.18) 

JdV 

It is worth noting that the resulting series features terms containing two or more consecutive factors of G N . As 
G N is diagonal, such terms correspond to electron-to-electron and hole-to-hole reflection processes. At first sight, 
the presence of such terms might be disconcerting, given the charge-interconverting character of Andreev reflection. 
However, it should be recalled that not only does the superconducting surround interconvert electrons and holes, but 
also it confines these quasiparticles to the normal region. For example, consider the series of terms that contain no 
superconducting Green functions G s . This series is precisely the Dirichlet series obtained in BB-I. Furthermore, this 
series can be embedded in any term of the MSE, which amounts to replacing the fundamental (i.e. unconfined) Green 
function by a suitably confined Green function. Therefore, terms involving consecutive factors of G N , rather than 
being disconcerting, are necessary to correct the free-propagation term, doing so by cancelling the Feynman paths 
that venture into the superconductor. 



C. Effective boundary conditions 

Before proceeding with our main issues (viz. the computation of the BDG Green function inside the billiard), we 
pause to pose and answer two questions: (i) Is there any boundary condition that can be imposed on G 11 so that G 11 
can be computed by solving the BDG Green function equation solely in V, i.e., without any reference to the region V. 
And if so, (ii) what is the precise form of this boundary condition? (If such an approach turns out to possible then 
one could dispense with the cumbersome task of dealing with Green functions having arguments outside V, as well as 
the concomitant need to match Green functions across the boundary.) 

To see that such a boundary condition does inde ed e xist, and to determine its explicit form, we substitute into 



Eq. (3.17) the parametrization of G 11 given in Eq. (3.1) in terms of fi n . Then all reference to fi n cancels, and we 



arrive at a nonlocal and billiard-shape-dependent effective boundary condition obeyed by G n , viz., 

G" = -2G S K°° cr 3 SG". (3.19) 

The reason that the boundary condition is nonlocal is that there exists the possibility of virtual propagation within 
the superconductor surrounding the billiard. The reason that the boundary condition is shape-dependent is that this 
virtual propagation outside the billiard is modified (from the value it would have in a homogeneous superconductor) 
due to the presence of the normal region. 
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IV. ASYMPTOTICS OF THE MULTIPLE SCATTERING EXPANSION 



Up to the present point, our investigation of the Green function for the BDG wave equation has been exact, and 
we have developed the machinery — the MSE — for computing this Green function in terms of the fundamental (N and 
S) Green functions and the shape of the billiard. The construction, however, is in terms of an infinite series, each 
term in which involves repeated integration over the boundary of the billiard. Thus, the direct computation of an 
arbitrary term in this series is prohibitively difficult, unless the shape of the billiard is exceptionally simple. To make 
progress we therefore need to invoke some approximation scheme and, as the Fermi wavelength is taken to be much 
smaller than the characteristic linear dimension of the billiard L, a very natural one to consider is the semiclassical 
approximation. In the present setting, this involves the evaluating of the repeated boundary integrals via a short-wave 
asymptotic approximation scheme. What we mean by this is that we seek an asymptotic approximation for every 
term in the MSE, the expansion parameter being 1/kpL, where kp is the Fermi wave vector; having invoked such an 
approximation, we shall re-sum the MSE. 

By following this scheme we shall be able to obtain, inter alia, the Green function for single-particle excitations, 
as well as a trace formula for the oscillatory part of the density of energy eigenvalues, valid in the semiclassical 
regime. What we mean by a trace formula is an explicit formula for the oscillatory part of the density of energy 
eigenvalues, expressed in terms of a sum over all closed semiclassical particle orbits. As we shall see, by virtue of the 
retro-reflective character of Andreev reflection, this sum over particle orbits is quite distinct from that arising in the 
setting of conventional billiards. 

Our approach has the virtue of delivering results not only for the DOS at the coarsest of energy resolutions (i.e. the 
Andreev level of approximation, in which motion is confined to chords traversing the billiard) but also at the finer 
level, thus revealing the mesoscale oscillations due to the quantal particle motion transverse to each chord. 

A. Classical dynamics in Andreev billiards 

The purpose of this subsection is to make a brief intermezzo in which discuss the physics of Andreev reflection from 
the point of view of geometrical optics. To this end, we develop stationary phase arguments aimed at elucidating the 
origin of the retro-reflective character of Andreev reflection. Along the way, we shall see that owing to the difference in 
the wavelengths of the incoming and reflected quasiparticles there is imperfectness in the this retro- reflection, i.e., the 
reflected excitation does not, in general, precisely retrace the path of the incoming one. Arguments of this type will be 
useful, subsequently, when we come to incorporate quantum fluctuations around the classical trajectories associated 
with Andreev reflection. 





h + 
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FIG. 3. Geometry for Andreev reflection from a planar surface. Initial and final points (i.e. the off-boundary points) are 
kept fixed, and the point of reflection from the boundary is then determined via the stationarity of the corresponding phase. 

It is well known that Eq. ( |2.5| ) gives rise to the Andreev reflection phenomenon, in which the electrons arriving from 
the normal metal are converted into holes (and vice versa) at the superconductor boundary, and are retro-reflected 
(i.e. have the excitation velocity reversed). In the present section we discuss the classical limit of this reflection process 
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by making use of the method of stationary phase (i.e. via the principle of least action). Throughout this section we 
make the (physical-optics-type) assumption that an electron wave having energy E and traveling a distance r acquires 
a phase 



whereas a hole traveling the same direction acquires the phase 



-ik— r 



(4.1) 



(4.2) 



where k + and fc_ are, respectively, the wavevectors appropriate for for particle and hole motion in the normal region. 
The energy dependence of these wavevectors is given by 



fc± = ^/fi±E. 



(4.3) 



Following t he s tandar d op tics-type approach, we envision some process (for an example see Fig. g) and then, by 
using Eqs. (hi) and (|4.2| ), we calculate the total phase acquired during this process. In the classical limit, the 
dominant process is the one (or ones) that make stationary this total phase, and hence determines information such 
as relationships between angles of incidence and reflection. Thus, in effect we are finding the rules of classical dynamics 
for Andreev billiards. 

At this point we have to introduce the physics of Andreev reflection "by hand,' and do so by requiring that after 
one reflection from the billiard boundary an electron is converted into a hole (and vice versa) . Thus we are demanding 
that there is no electron-to-electron (or hole-to-hole) scattering (due to a single reflection). Under these conditions, 
any scattering process can be analyzed in terms of the basic electron-to-hole and hole-to-electron processes. We need 
only examine one of these because the corresponding phases are identical and thus, at stationarity the two processes 
have the same geometry. In other words, the stationary path describing an incoming electron and scattered hole may 
be reversed (by reversing the direction of propagation of each particle) to give the stationary path of an incoming 
hole and the scattered electron. 




FIG. 4. The analogy between Andreev reflection and optical refraction 



It might be useful to note the similarity between the phenomena of Andreev reflection and the the refraction of 
light. The common feature is that, in both settings, before and after scattering the waves have the same frequency 
but differing wavevector magnitudes. In the case of the refraction of light, the wavevector is changed because the 
wave enters a medium with a distinct index of refraction. In the case of Andreev reflection, the wavevector is instead 
changed because the electron wave is converted into a hole wave, the latter having a distinct dispersion relation. In 
fact, by reversing the sign of the phase for a hole wave, as well as the direction of propagation of the hole wave, one 
transforms the Andreev reflection process into the familiar optical refraction process. Thus, one has an electron-to- 
hole reflection law that is essentially identical to Snell refraction (rather than specular reflection). Thus, Andreev 
reflection looks like optical refraction but with the outgoing direction reversed with respect to the scattering point 
(see Fig. |). 
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To quantify these remarks, consider the problem depicted in Fig. ||, in which an electron arriving from a fixed point 
(1) is reflected and converted at the variable point x into a hole, which then propagates to another fixed point (2). 
The classical path corresponds to the value of x at which the total phase for the process is stationary with respect to 
variations of x. For the process at hand, total phase is given by 



k+Jhl + x 2 - k^Jh 2 _ + {x-l) 2 , 



for which the stationarity condition reads 



x — I 



yjh\ + X 2 ^]h 2 _ + {X- l)> 



= 0. 



(4.4) 



(4.5) 



By rewriting this condition in terms of the angles of incidence and reflection (i.e. 9 + and 
recovers the Snell's law form: 

fcj. sin 6i = fc_ sin . 



shown in Fig. one 



(4.6) 



By using Eq. (4.6) one can construct the stationary paths for a n A ndreev billiard, just as one does in the case of 
geometrical optics. W hen there is more than one reflection, Eq. (4.6) must be satisfied at each one. 

One feature of Eq. (4.6) is that it makes evident the fact that the reflected particle is not, in general, perfectly 
retro-reflected (i.e. 9+^9-). However, when k + and fc_ are very close to each other, 9+ and 9- will be, too. Let us 
now calculate this small deflection angle 9 - — 9 + in terms of \x and E. To do this, let us assume that 9+ and 6L are 
close and that E/fi-^i 1, and expand Eq. ( |4.6| ) to obtain 



(E/n) tan6» 4 



(4.7) 



As expected, the deflection angle is 0(E/n), unless the incident direction grazes the boundary. This qualification 
divides the space of incoming trajectories into two classes: (i) a large fraction, occupying most of the phase spa ce, i n 
which tan#+ is of order unity, and (ii) the rest, in which the deflection angle 9- — 9 + is not small. Equation ([4.7]), 
although approximate, provides a guide for addressing whether or not deflections (i.e. imperfectness in retro-reflection) 
needs to be taken into account. 




FIG. 5. Typical closed classical trajectory in an Andreev billiard, 
paths. The imperfectness of the retro-reflection is exaggerated. 



Black lines depict electron paths; gray lines depict hole 



We are now at the p oint where we can construct classical dynamics in an Andreev billiard. For any trajectory, the 
reflection rule in Eq. (4.6) has to be satisfied at every reflection point. This will generate a type of dynamics that 
differs from that generated by the usual specular reflection rule (in which the outgoing angle equals the incoming 
one). An example of a closed trajectory in an Andreev billiard is shown in Fig. |[ 
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B. Asymptotics of the fundamental Green functions and their derivatives 



In the present section we investigate the asymptotic behavior of the fundamental Green functions G N and G s for 
both small and large values of their (position) arguments. This investigation will allow us to estimate the relative 
dominance of various processes, and thus to organize the multiple scattering expansion for the exact BDG Green 
function into a form suitable for establishing its approximate behavior at large kpL. 

The asymptotic behavior of G N ' S is related to the corresponding Helmholtz Green function g± , which can be 
represented through the Fourier integral 

N.s m f d d p e'P' 

■ 9± = v^ 1 ' (48) 

N S ... • • N S 

where k = k± , depending on the Green function in question. Then the asymptotic behavior of g± for large kl can 
be obtained from the asymptotic evaluation of this Fourier integral, which gives 

N,s m , ■ ( *±' 3 \ ^ exp(± ? fc^ s Z T *7r(rf-l)/4) Q . 
fc'(0«±*^J ^ ■ (4-9) 

The derivatives of the Green functions for large kl can be obtained by differentiating this asymptotic expression. 
Determining the small kl asymptotics of g±' S (l) is more tricky. By scaling p with I we get 

•'<<>= ,<2 -VIf^w- (4 ' 10) 

N S 

where a = pl 7 and g is shorthand for any of the four g± . Thus for small kl we have 

(1) „ f if d^2; (4n) 

gy) \ln(M), iid = 2. [ ' 




FIG. 6. Approximating the surface by its tangent plane 

Having determined the form of g for small kl, we now investigate what can be said about g{l), dg(l) and dSg(l) 
when 1 is a vector connecting two nearby points on the surface dV. The geometry (for the d = 3 case ) is illustrated 
m Fig. |. In this fi gure, 1 = j3 — a, where (3 and a are points on the surface. In the following, we shall work 
in d = 3, although results obtained will be applicable for all d. Without loss of generality, we choose a to be the 
origin of our coordinate system, with the z direction coinciding with the inward normal direction of the surface at 
a. The remaining directions are chosen arbitrarily (at least for the time being), the only constraint being that the 
coordinate-system be right handed. In this coordinate system the surface may be defined locally through an equation 
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of the form z = f(x, y), where x and y span the tangent plane, and dxf\^ x y )— (o o) = ®nf\(x y)—(o o) = ^' Then a P om t 
(3 (near a) on the surface will approximately have the coordinates 

However, by a suitable rotation within the tangent plane the equation for f3 may be written as 

where R\ and i?2 are the two principal radii of curvature of the surface at the point a.. For reasons that should soon 
become clear, we define the vector (3// to be the projection of (3 on to the tangent plane at a. Next, we make two 
assumptions: 

1. the radii of curvature are on the order of the linear size of the billiard, i.e., Ri^/L = 0(1). 

2. the region we of interest around a has a linear size on the order of kp 1 . 
Under these assumptions, we have that 

l=\(3-cc\ = \{3 // \{l + 0{(k F L)- 2 )}. (4.14) 

Let us now turn to dg. In the coordinate system specified above, the normal vector ng at the surface point (3 is 
given by 

n ^ = (^'^' 1 ) + ° ((fcFLr2) - (4 ' 15) 

Then the quantity dl is given by 

{{kpL) ^ = ° (i) ■ (4 - i6) 

Hence, we have that the normal derivative of the Green function is given by 



dl I \Ri R 2 
By generalizing to arbitrary dimensionality d we obtain 

dg{l) = O (L- x l {2 -^y (4.18) 

In order to evaluate d5g, we consider a third point on dV, which we denote by 7, focus on the quantity 381' (where 
V = |/3 — 7I) and, at the end of our calculation, let 7 tend to a . We shall use primed coordinates x' and y' for 7. In 
this way, we find that 

°« ■ <■* • ™» • ^ - * - (is s) - 1 

32 



= (^7 V(* - ^') 2 + (y - y') 2 + (* - ^') : 



„2 „/2 „,2 „J1 \ 2 1 1 / / 7/2 



i^-^s^y . . (4 , 9) 



z' 3 v. 2i?i 2i? 2 y z' /'V V^ 2 

We are now in a position to evaluate dSg, for which we find 

ai9(i) .| M + || (ai)TO ._I|( 1+0 (;)). (4 , 0) 

In the MSE, for each term that includes the product g ddg there is a corresponding term in which g dSg is replaced 
by dgdg. It is therefore desirable to estimate relative size of these terms for small values of their arguments. The 
asymptotic formulas given in the prcssent section are useful for this comparison, giving 

l*(0*ffll-W)fl«.(0|x{2(^ w)i (4.21) 
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C. Asymptotic expansion for the quantal amplitude 



The MSE of Sec. Ill has provided us with a series expansion for the Green function (i.e. the full quantum-mechanical 
amplitude for the propagation of a quasiparticle excitation from one point x to another x'). This series expresses the 
correction to the free-space propagation of quasiparticles caused by multiple scattering processes of the quasip articles 
from the boundary that separates the normal and superconducting regions of the billiard. A generic term features the 
following possible processes: (i) inner reflections (which are marked by the occurrence of an adjacent pair of normal- 
state Green functions in the algebraic expression for the contribution), (ii) transmissions (marked by the occurrence 
of an adjacent pair of Green functions, one normal and one superconducting), (iii) and outer reflections (marked by 
the occurrence of an adjacent pair of superconducting-state Green functions). Throughout the present section, wc 
shall assume that d — 3. (The extension of the following discussion to general d is straightforward.) 

The generic contribution to the amplitude involving a total of n reflections and scatterings can be written as 

„4(x,x') = J M{x,ai,a 2 , ■ ■ ■ , a„,x') expik F S(x,ai,a 2 , ■ ■ ■ ,a„,x'), (4.22) 
where the integral is taken over all values of {c*i, ct2, ■ ■ ■ , ot n -\, c* n }, each element ranging over the surfac e dV. Th e 



modulus function A4(x,c*i, 012,1x3, ... , a n , x') is real and, as can be seen from the iterative solution of Eqs. (3.1 -3.3d) 



is a sum of products of functions such as \oti — a?i+i| , first and second normal derivatives of this function, as well 
as a polynomial in kp and /c^' S . The phase function S(x,oti 7 ot2,ot 3 , ■ ■ ■ ,ot n ,x') is, in general, complex and, as can 



also be seen from the iterative solution of Eqs. (3. 1-3. 3d), is a sum of terms each of the form (fc^' /A;f) \oti — c*i+i|. 

The approximation scheme that we shall invoke involves the short-wave asymptotic expansion of this quantal 
amplitude A, valid for large kpL. The method used for the construction of this expansion is the asymptotic expansion 
of the multiple integrals appearing in the terms in the MSE. (See, e.g., Ref. [g9| for a discussion of the asymptotic 
expansion of multiple Fourier integrals.) This method allows one to approximate -4(x, x') as a sum over critical points 
(c.p.) of the domain of integration (which, in this case, is the 2n-dimensional manifold 971 = dV x dV x • • ■ x dV): 

.4(x,x')«^A. P .(x,x'), (4.23) 

c.p. 

where „4 c . P .(x, x') depends solely on the local properties of M. and S at the critical point. As kpL — •> 00, corrections 
to this formula vanish faster than any power of (kpL)^ 1 |pC}| . In other words, for large kpL contributions from the 
critical points dominate the total amplitude. A critical point {pt\, a 2 , ot^, • • • , ct^) can be interpreted as a trajectory 
(although not necessarily a classical one) in which an excitation travels from point x to x', along the way scattering 
at the points a\, a^, ■ . . , c*^, so that A c . p .(x, x') can be interpreted as the amplitude corresponding to this trajectory. 

Before actually proceeding with the construction of the expansion of the multiple integrals appearing in the terms 
in the MSE, we first classify the critical points of 971: these are 

1. Points at which the gradient of the phase, i.e., V ai 5(x, c*i, a 2 , a 3 , • • • , a„, x') vanishes for all i. 

2. Points at which M.{x, ol\, a 2 , • ■ ■ , ot n _x, a n , x') has a singularity. 

3. Points at which M. or S fail to be infinitely differentiablc. 

4. All points on the boundary of the manifold 97t. 

5. Points satisfying criteria (1-4) in a mixed sense, i.e., points satisfying criterion (1) within a submanifold of points 
satisfying criterion (2) within a submanifold of points satisfying criterion (3) within a submanifold of points 
satisfying criterion (4). 

In the present setting, 971 has no boundaries and, thus, there are no Type 4 critical points. As for Type 3 critical points, 
when dV is infinitely differentiable, so are A4 and 5 and, hence, there are no Type 3 critical points either Thus, 
the only possible types of critical point are (1), (2) and (5). In present case, M. consists of products of functions such 
as \a.i — and its first and second normal derivatives. Thus, critical points of types 2 and 5 occur whenever 

one or more of the propagation distances \ati — £* 2 :+i| vanish. In accordance with the trajectory interpretation of 
critical points, in which the sequence {c*i, 012, • • • , ot n } defines the trajectory, we call the part of the critical trajectory 
having vanishing propagation distance a zero-length path. Then, all Type 5 critical points can be generated from 
Type 1 critical points by the insertion of zero-length paths. Stated technically, if (oti, a 2 , • • • , c*i, CKi+i, • • ■ , c* n -i) is 
a critical point of the [(2n — 2)-dimensional] manifold 371' then the point (qi, O62, ■ ■ ■ oti, oti, cci+i, • • • , ot n —i) will be a 
critical point of Type 5 in the [2n-dimensional] manifold 971. Criterion (1) amounts to the familiar stationary-phase 
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approximation for the amplitude function, because the points at which all gradients of S vanish are the stationary 
phase points. 

Throughout this Paper we are interested in leading- order contributions to quantal amplitudes. Thus, it is useful to 
determine whether or not a process contributes to the full amplitude at leading order. To do this, we must be able 
to estimate the order-of- magnitude of contributions from different types of critical point. To this end, let us consider 
a type 1 critical point (i.e. a trajectory defined via the principle of stationary phase), and a Type 5 critical point 
constructed from this Type 1 critical point via the insertion of a zero-length path. For the sake of simplicity, let us 
consider as our Type 1 critical point a very simple amplitude, i.e., one having just one reflection: 

A(x, x') = j dg(x, a) g(a, x'), (4.24) 

where g is a generic Helmholtz Green function. [For the purposes of determining the order-of-magnitude of the 
contribution from various critical points, whether the Green function is g N or g s is irrelevant.] By using the asymptotic 
formulas for g presented in Sec. IV B, it is possible to write asymptotically (up to from numerical factors) 

.A(x,x')~ [ -. £ -exp(ik\x-a\+ik'\a-x'\). (4.25) 

J |x — a\ \a — x'\ v ' 

Here, k and k' can be ±fc+'_. The position (a c ) of the critical point will depend on the chosen values of k and k' via 
the stationarity condition. However, owing to the facts that k and k' are both O(kp) and that |x — a c \ and \a c — x'| 
are both O(L), it is possible to estimate the order-of-magnitude of the contribution associated with (a c ) to be 

„4(x,x') = O (k F /L 2 ) exp (ik\x- a c | + ik'\a c - x'Q J dxdy exp (ik F L(Ax 2 + Bxy + Cy 2 )), (4.26) 

where we have expanded the phase to second order in deviations from (a c ). From dimensional considerations we 
know that A,B and C are all (D(L~ 2 ) and, thus, that 

A{x, x') = O (1/L) exp (ik\x - a c \ + ik'\a c - x'|) . (4.27) 



Now let us consider the insertion of a zero- length path. The rules described in Sec. 1IIB for constructing a term in 
the MSE allow three possible types of such insertions: (i) the insertion of g, (ii) the insertion of d6g, and (iii) the 
insertion of dg. More specifically, we are interested in the amplitudes 

A(x, x') = J dg(x, a) g(a, (3) 6g((3, x'), (4.28a) 

An{x,x') = J g{x,ct)d5g{a,/3)g((3,x'), (4.28b) 

Au(x, x') = J dg{x, a) dg{a, (3) g(J3, x'), (4.28c) 

and their leading-order asymptotic contribution due to the critical point (oc,f3) = (a c ,ot c ). Let us start with Ai. By 
using a coordinate system centered at a c and the short-distance asymptotics of g, and expanding S to second order 
around a r , we find that 



Ai(x,x') = O ( % ] e 4 ^ / da a dan ^LjAxl+B^yMl+A'^+B'^y.+Cy^ 

\L 2 J J P J(x a - x d Y + (y a - y 8 Y 



yj{x a - xp) 2 + (y a - yp) 



O (^) O (J^j e lS < J dxdy exp (ik F L(A"x 2 + B"xy + C"y 2 )) = O (~\ 



e lb % (4.29) 

where A, A', A", B, B', B", C, C and C" are 0{l/L 2 ), D and E are O(l), and S c = k\x - a c | + k'\oc c - x'|. Thus, 
Ai contributes at the same order as A. A similar calculation provides an estimate of the order-of-magnitude of An, 
and show it to be of the same order as A and A\. Next, let us consider By using the short-distance behavior 

dg{l)*o(\\^lKo(\\g{l), (4.30) 



L) dl \L r 
the order of magnitude of Am can be deduced from the estimate of A 
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(4.31) 



Thus, we see that Am contributes to the full amplitude only at subleading order. 

In the usual case of a Schrodinger billiard with hard walls (i.e. Dirichlet boundary conditions), the asymptotic 
contribution coming from the critical points of Type 2 and Type 5 have been shown to be smaller than the stationary- 
phase (i.e. Type 1 critical points) contributions, and by a factor of k?L ||. The reason for this is that in the hard 
wall case all the singularities of M. are due to dg's. Thus, for the Schrodinger billiard with hard walls, as far as the 
leading-order contribution is concerned the Type 2 and Type 5 critical points are irrelevant. Thus, in such billiards, 
the leading asymptotic contribution comes from the stationary-phase points. 

In constrast with the case of Schrodinger billiards with hard walls, in Andrccv billiards (and Schrodinger billiards 
with soft walls, i.e., with a finite bounding potential) the Type 2 and Type 5 critical points do not necessarily give 
only subleading-order contributions. More specifically, in Andreev billiards, in addition to its dg singularities, M can 
have additional singularities due to g and dSg. As shown above, both of these singularities contribute at leading order 
in the (fcpL) -1 expansion p^ |. Thus, Type 2 and Type 5 critical points are relevant for Andreev billiards. 

Having determined the significance of Type 5 critical points for Andreev billiards, we now examine these critical 
points more closely. The first important observation is that the insertion of a zero-length path does not change the 
value of the phase function S. As all Type 5 critical points can be regarded as originating from Type 1 critical points 
via the insertion of a suitable number of zero-length paths, the phase of any Type 5 critical point is equal to that of 
the originating Type 1 critical point. The second important observation is that, because the phase is not changed by 
the insertion of zero-length paths, all amplitudes originating (via insertions) from a given Type 1 critical point carry 
a common phase, and therefore add coherently. Thus, the effect of the Type 5 critical points is to modify (but not 
necessarily increase) the amplitude of the contribution to the originating Type 1 critical point. 

In order to make less abstract the issue discussed in the previous paragraph, consider the example of reflection 
from an infinite plane boundary in the short-wave asymptotic limit. For the case of the hard Schrodinger billiard 
there is a single critical point, which is of the stationary-phase type: it is the classical reflection point (for which the 
angle of incidence equals the angle of reflection). For the case of the Andreev billiard there are two possible electron 
reflections: electron-to-hole and electron-to-electron. These two processes have differing phases and, correspondingly, 
differing stationary-phase (i.e. classical reflection) points. If one were to take into account only the stationary-phase 
(i.e. Type 1) points then one would find that the amplitude for electron-to-hole reflection would vanish, whereas that 
for clcctron-to-clcctron reflection would be of order unity. However, this finding would be misleading, owing to the 
fact that the set of critical points that contribute at leading order (in the the short-wave asymptotic limit) is much 
larger. To see this, focus on the case of electron-to-electron reflection. Let us label the classical reflection point by ot c . 
Then set of critical points is (a c ), (c*c, c*c), (c*c, c*c, <*c), (c*c, &c, c*c, £*c), etc., i.e., there is the possibility of multiple 
scatterings from the boundary, all taking place in the vicinity of the classical reflection point. These additional 
critical points correct the amplitude for this scattering process, and yield the expected result, namely that the net 
electron-to-electron amplitude is very small. The origin of this correction, then, is that multiple virtual propagations, 
inside the superconductor but near to the classical reflection point, decrease the amplitude of the electron-to-electron 
reflection process. Mutatis mutandis, this mechanism increases the electron-to-hole reflection amplitude, leading to 
the familiar physics of Andreev reflection. 

There is a simple physical explanation for this mechanism. Quasiparticle propagation in the bulk of the normal 
region is not affected by the superconducting surround, except via those Feynman paths that pass nearby the boundary. 
In the short-wave asymptotic limit, this occurs near reflection points. Thus, quantum mechanically, there is an effective 
volume around the boundary in which propagation is modified due to the amplitude for electron-to-hole conversion 
(and vice versa). Hence, there is an effective volume around the classical reflection points, and in this volume 
multiple scatterings convert electrons arriving from the interior of the normal region into holes departing for the 
interior. Classically, the volume for such processes is zero, i.e., the conversion takes place precisely at the reflection 
point. Thus, zero- length propagation at the boundary is responsible for the electron-hole interconversion aspect of 
Andreev reflection, whereas the requirement of phase-stationarity, applied to propagation in the interior of the normal 
region, is responsible for the retro-reflection aspect of Andreev reflection. 

D. Integrating out propagation in the superconducting region: Effective reflection 

We now actually evaluate the contribution to the short-wave asymptotic approximation to the Green function that 
arises from all critical points involving zero-length propagation. In doing this, we collect contributions from Type 5 
critical points, and arrive at the expected result that reflection leads to almost complete electron-hole interconversion. 
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We start with the expression (3.17) for /x" which, for the sake of convenience, we rewrite here along with the explicit 



form of K°° obtained from definition ( 3.14b ) 



H [i = -20-3 G N + 4o- 3 G s K°° cr 3 <5G N + {-2<r 3 <9G N + 4er 3 G s K°° cr 3 <9<5G N } /x", 
K°°= (l + 2er 3 .<5G s ) _1 . 



As we have shown in Sec. IV C, in the short-wave asymptotic limit, critical trajectories that have been obtained from 
a Type 1 critical trajectory (i.e. a pure stationary-phase trajectory) by the insertion of zero-length propagations of 
dSG N and G s contribute at the same order as the original Type 1 critical trajectory. Thus, such contributions should 
be summed to all orders. On the other hand, critical trajectories obtained from Type 1 critical trajectories by the 
insertion of <9G N and 5G S contribute only at sub-leading orders. Thus, it is appropriate to ignore such contributions. 

Moreover, we are considering situations in which the range of G s is much smaller than the size of the billiard. 
Thus, all critical points that include finite-range superconducting propagation are suppressed exponentially (in the 
size of the billiard), despite their being formally of leading order. We shall therefore neglect them, at least for the 
time being. Such contributions constitute the single-particle tunneling amplitude through the classically-inaccessible 
S region. Below, in App. |c], we shall study the consequences of relaxing the condition that the range of G s be much 
smaller than the size of the billiard. We shall then show that in settings involving convex billiards (i.e. billiards for 
which all chords lie inside the billiard) such contributions cancel each other at leading asymptotic order. Thus, for 
the purposes of a leading-order calculation we can make the approximation <5G S 0. It follows that 

K°° « I. (4.32) 

The only remaining appearances of the superconducting Green function (i.e. G s with no derivatives) in the MSE 
generate critical points having zero-length superconducting propagation. Moreover, the only Green function that 
contributes at leading order to both zero-length and nonzero-length propagation is <9<5G N . In order to re-sum all 
possible zero-length propagations it is natural to separate the operator <9<5G N into two parts: one solely generating 
zero-length propagation; the other solely generating finite- length propagation Q . A convenient way to do this is to 
define the following operators: 

(d6G™ F) (a) = I da d a dpG N (cx,l3)F(l3)w(a - 0), (4.33a) 
JdV 

(dSGf¥){a)= [ da d a d,3G N (a,/3)F(/3)(l-w(cx-/3)), (4.33b) 

JdV 

where w(a — 0) is a smooth function that equals unity whenever a and (3 are close to one another and vanishes 
whenever a and (3 are far away from one another. The effect of this function is to isolate the critical point at (3 = a 
from the remaining critical points, the latter having finite-length propagation involving cWG N . For the purposes of 
our asymptotic expansion, the isolating function w is only a convenience, and its particular form does not affect the 
final results (as long as the range of w precludes its enveloping simultaneously any pairs of critical points). Then the 
equation for /x 11 becomes 

/x u « -20-3 G N + 4<t 3 G s cr 3 <5G N + {40-3 G s 0-3 dSG™ + (-2er 3 <9G N + 4<r 3 G s <x 3 dSGf) } /x". (4.34) 

Having decomposed the kernel of this equation into two pieces (the first consisting of critical points having zero-length 
propagation and the second consisting of critical points having finite- length propagation) we invert this equation with 
respect to the former piece, obtaining 

// sa (I - 4<r 3 G s cr 3 ^G^) _1 {-2<t 3 G N + 4<r 3 G s <x 3 <SG n + (-2<r 3 <9G N + 4<r 3 G s <x 3 dSGf) /x"} . (4.35) 

We now define the renormalized Green function G R (a,x'): 

G R (a,x') = (l-4cr 3 G s cr3 5(5Gf) _1 { -<x 3 G N + 20-3 G s <t 3 <5G n } . (4.36) 



In terms of G R , Eq. ( 4.35 ) becomes 

fi li « 2G R + 29G R /x". (4.37) 

We note that this equation, no critical points containing zero-length propagation contribute to /x 11 at leading order. 
Thus, the summation of short-range critical orbits is achieved via the calculation of G R in the short-wave asymptotic 
approximation. 
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The main contributions to the integrals implied in Eq. (4.36) come from the neighborhood of a. This follows from: 
(i) the fact that isolating function is short-ranged; (i) the fact that G s is finite-ranged; and (iii) the assumption that 
the billiard is large enough to exponentially suppress any finite-range critical points produced by G s . Therefore, it is 
adequate to approximate the boundary surface dV around a. The lowest-order approximation to dV is the tangent 
plane at a. The corrections to this approximation are smaller by a factor of (&F-R) -1 , where R is the radius of 
curvature at the point ex.. Throughout this Paper we are assuming that the surface dV is sufficiently smooth that R is 
of order L, i.e., the radius of curvature is comparable to the billiard size. Thus, corrections due to the curvature of the 
surface do not contribute at leading order. Having replaced dV by a tangent plane, the integral equation for G R (a, x') 
becomes solvable, owing to the resulting translational invariance in all directions parallel to the tangent plane. Thus, 
by introducing the two-dimensional Fourier transform (2DFT; see App. 0) of all Green functions appearing in the 



integral equation (4.36), it is straightforward to obtain the following algebraic result for the 2DFT of the renormalized 
Green function: 



G R (p, z') = {I — 4<r 3 G s (p) <x 3 d5G N (p)} 1 {-<x 3 G N (p, *') + 2<x 3 G s (p) <x 3 <5G N (p, z')} 



(4.38) 



to which there are corrections of order (fcpi?) . (This Grenn function is exactly the Green function for a planar 
boundary.) Here and elsewhere, p denotes the magnitude of the 2D vector p conjugate to the position-vector in the 
tangent plane. 

We now embark on the task of inverting the 2DFT G R (p, z') in order to obtain the approximate real-space renor- 
malized Green function G R (c*,x). Thus, we need to evaluate the integral 

G R '(a, x') = J J^L G R (p, z') exp ip ■ (a - /3), (4.39) 

where (3 is the component of x parallel to the plane and z' is the perpendicular component. The quantities G N (p, z') 
and SG N (p, z'), needed to construct G R (p, z'), are derived in App. [B], where they are found to be given by 





SG (p,z) = ( n i ) = (n ill o e -a-{p)\z'\ J, ( 4 -40b) 



in which a±(p) = ^Jp 2 — k±, the square roots being taken such that their real parts are always positive. Note that 

the only z' dependence in G R (p, z 1 ) comes from G N (p, z 1 ) an d i5G N (p, z'), and it is only in these terms that there is 
exponential dependence on p. By inserting Eqs. ( 4.40a ) and ( 4.40b| ) into Eq. ( 4.3§| ) we find that G R (a,x') is given 



by 

f j2 / I c - a +(p)l 2 'l \ 

G R (",xO = y^R(p)f ^ -^L_e--W^lJ eXpiP ' (a " /3) ' 
Here, R(p) is a certain algebraic function of p, and is defined by 

R(p)^{l-4 ( r 3 G s (p) ( r 3 dSG N (p)y 1 j(j + 2a 3 G s (p) _ a ° (p) )}> (4.41b) 

The analytic expression for R(p) is obtained using the 2DFTs of G s (p) and dSG N (p): 

s 1 f E A 1 / 1 1 \ 1 / 1 1 

[P) ~ 2\i^A*-E2 zVA^-E^ 2 } \2a s + (p) 2a s (p)) + 2 °" 3 \2a*(p) + 2a s _(p) 

dSG N (p) = - i I (o+ (p) - a_ (p)) - i a 3 (a+ (p) + a_ (p)) . (4.42b) 

We evaluate the integral in Eq. ( |4.41a| ) via the method of stationary phase, which becomes exact in the limit kp\a — 
x'| — > oo. From the form of a±(p) we see that for values of p for which a±(p) is essentially imaginary (i.e. for 
p < Reefc-t) there exists the possibility of a stationary-phase point. Note, however, that a + behaves differently from 
a_, due to the fact that the imaginary parts of k\ have opposing signs. For p < Refc_ (note that k- is always smaller 
than k + ) we have 
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a±(p) = =R\/ fc ± -p 2 . 
The stationary-phase point is defined by the condition 

from which we see that stationary-phase point p c satisfies 

(a - 0) , p, 



±- 



k± - |p c 



(4.43) 



(4.44) 



(4.45) 



and that |p c | 2 has the value k± sin 2 9 ax ' , where 9 ax ' is defined to be the angle between the normal vector at the 
surface point a and the vector x' — ex.. Then the effective Green function can be asymptotically approximated as 



G R (a,x') 



R++ (k% sin 2 e ax . ) R+- (k 2 _ sin 2 
R-+(k 2 + sm 2 9 ax ,) #__(fc 2 sm 2 



d 2 p 



} -a + (p)\*'\ 

2a+(p) 





3 ip-(a-/3) 



R++ 



5+ (a — x') 

-#_(a-x') 



(4.46) 



The second line is obtained by noting that the integral is, in fact, the 2DFT of G N (see App. |b|). The amplitudes R++, 
R-i , R |_ and R can now be respectively interpreted as t he elec tron-electron, electron-hole, hole-electron and hole- 
hole reflection amplitudes, and can be obtained from Eq. (|4.41b ). These amplitudes are, in general, nonvanishing. 

However the charge-preserving amplitudes (i.e. R++ and R ) are smaller than the charge interconverting amplitudes 

(i.e. R-\ and R (.) by a factor of A//xcos 2 6>. In order to evaluate R to leading order in A/^tcos 2 #, we first note 

that for 9 ax i not near 7r/2 (i.e. not near grazing) one has the following approximations for a(p): 



a±(k~F svat 
a±(k-p smt 



^fik-p cos( 
^fik-p cos( 



O(A /m) 



By applying these approximations to Eqs. (4.42a-4.42b) we find 

( E A 



G b (fc F sin( 



dSG N {k F sin6» Q2 



f 



: <T 2 



(4.47a) 
(4.47b) 



(4.48a) 



By using these expressions in Eq. (4.41b) we obtain R(fcp sm9 ax i) 



gl(a,x' 



2«a/Aq - E 2 2i^J Aq - E 2 'J -ik F cos 9 a 

(4.48b) 

e~ l<p (Ti and, thus, from Eq. (4.46Q we obtain 



-lik-p cos 9 a 



G R (a,x')«cxp(-^ + i|) ( 



9-(<x, x') 




(4.49) 



The off-diagonal structure represents the total electron-hole interconversion that occurs for large perpendicular mo- 
menta. However, strictly speaking, electron-hole interconversion is not perfect, i.e., the amplitudes R++ and R do 

not vanish. Moreover, these charge-preserving amplitudes increase, as the angle of incid ence appr oach es tt/2. Note 
the difference between the regimes of validity for the approximate expressions Eq. (4.46) and Eq. (4.49): the former 
becomes valid for k-pL ^> 1 whereas the latter becomes valid for k-pL ^> 1 and A/p, <C 1. However, in cither case the 
MSE can be cast into the following effective Multiple Reflection Expansion: 

G u (x,x') = G N (x-x') + 2 / da Q (9 ct G N (x~a)G R (a,x')+4 / da Q da <9 Q G N (x - a) d G R {a 7 0) G R (/3, x') 



i)V 



+8 



/ da a da dcr^ <9 Q G N (x - a) d p G K {a, (3) a 7 G R (/3, 7) G R ( 7 , x') + 

JdV 



(4.50) 



We stress two points about Eq. fl4.50| ): (i) it is free of leading-order short-range critical points [i.e. the contributions 
of short-range critical points to G u are smaller, by at least a factor of 0(1/ kpR), than the leading-order contribution, 
and this is what we were aiming for]; and (ii) all propagations inside the superconductor have been integrated out, 
leading to the effective reflection expansion for G 11 . 
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V. DENSITY OF STATES OSCILLATIONS 



In Sec. IV D we integrated out superconducting propagation by evaluating the short-range critical points in the 



MSE and, hence, we obtained an effective expansion for the Green function, which we have termed an MRE. In the 
present Section we shall focus on the density of states p(E), expressing this quantity in terms of the Green function 
which, in turn, w e expr ess via the MRE. In this section we shall ignore the effects of normal reflection, returning to 
them only in Sec. VB 4. Thus, by using the results of the previous section, we have 



p{E) 



1 




(x,x')+ 5 *(x,x') + 4 / dg"( X , a )dg R ( a ,(3)g*(l3,x') + dg ] 



l(x,a)dg*;(a,l3)gl 



(dgldgKdgKdgK---dg R g*:) 



[ (dgidgKdg R dgK---dg*g^ 

JdV 



(/3,x') 



(5.1) 



where terms with odd numbers of reflections van ish, owing to the off-diagonal structure of G R . First, let us note that 
the first two terms on the right hand side of Eq. (5.1), which have zero-length propagation and hence vanishing action, 
do not introduce any oscillations into the DOS. As these terms do not involve any surface integrals (and, hence, do 
not involve any surface effects), they produce the bulk DOS of a homogeneous N region Q: 



d d xlmg^{x - x 1 ) 



S 



d-l 



2(2tt) 



Vk 



d-2 



(5.2) 



where is the (d — l)-dimensional surface area of a <i-dimensional unit sphere and V is the volume of V. We now 
deal with the (remaining) critical points that contribute at leading order. These critical points are closed classical 
trajectories consisting of the propagation of quasiparticles through the bulk of the billiard (i.e. the N region), connected 
by reflections from the billiard boundary (i.e. the N-S interface). 

We shall distinguis h b etween two asymptotic approximation schemes for p, both of which are obtained by evaluating 
the integrals in Eq. (5.1) within the stationary-phase approximation, valid for for large k-pL and small A/ p. From a 
technical point of view, the difference between the two schemes concerns the stationary-phase points they use, which 
must be in accordance with the particular limits assumed for the parameters kpL and A/// (which the approximation 
becomes exact). 

Scheme A: The first scheme is, in essence, equivalent to the (by now conventional) adiabatic approximation to the 
wave function, first introduced by Andreev 0]; it becomes exact when energy-level spacing goes to zero, which occurs 
for the following limit: 



k-pL 



CO, 



A//x — > 0, and LA/kp — > constant. 



(5.3) 



In this scheme, an excitation undergoes perfect retro-reflection (i.e. perfect velocity-reversal) because in this limit 
the difference between k + and fc_ is ignored in the calculation of the critical trajectories. The resulting classical 
dynamics is confined to the chords of the billiard and, thus, is integrable, regardless of the shape of the billiard. 
However, for finite values of the parameters (i.e. large but finite kpL and small but finite A/ p), Scheme A produces 
the locally-energy-averaged DOS, which becomes numerically accurate only around the DOS singularities that it 
correctly captures. However, it fails to capture the DOS oscillations arising from the confinement of quasiparticles to 
the billiard. To capture these oscillations is the main motivation of the following scheme. 

Scheme B: In this scheme we shall take into account the previously neglected difference in electron and hole wavevec- 
tors. This leads to impcrfcctncss in retro-reflection because, upon reflection, the wavelengths of the incoming and 
outgoing waves are no longer identical (as happens with refraction except, of course, that the waves are now on the 
same side of the boundary). Technically speaking, approximation scheme B becomes exact when 



kpL — > oo, and A//i is a small parameter. 



(5.4) 



Now the classical trajectories are determined by the reflection rule given by Eq. (4.6) and, consequently, the dynamics 
is no longer a priori integrable; on the contrary, it is weakly chaotic for most billiard shapes p5[ |. 

In order to understand the distinction between asymptotic schemes A and B, consider the phase function S for a 
process having 2n reflections: 



S(x, ai,a 2 , • • • , ot 2r , 



k—£r, 



k-l 



C*2ti — 1 ,Q>2t) 



(5.5) 
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where £ a ^ ai+ i = |ctj — ctj+i|. Notice that S can be separated into two parts, a large one SAnd and a small one Si mp 
so that S = SaucI + Si mp , where 

SAnd(«^; C^2; ; ^2n) — 2 ^ai,a2 ~l~ ^-a2,c*3 ^3,(14 ' ^ct2?i — 1 ,c*2n ^ct2 T i,:c) 7 (5.6a) 

fc + — fc_ 

Simp(«£; t^2; i ^2n) — 2 ~t~ ^ai,a2 ~t~ ^a2,c*3 ~t~ ^3,04 ~t~ ^ct2?i — 1 ,c*2n ^ct2 T i,:c) • (5.6b) 

For n not too large (i.e. n <§; kpL), S'And = Oikph) and S; mp = O(LAfkp). Although it is clear that S'And is large 
and must be included in any stationary phase calculation, whether or not S lmp should be included in a stationary 
phase calculation depends precisely on the nature of the limiting scheme. If one solely uses SAnd for the calculation of 
stationary-phase points then the critical trajectories feature retro-reflection; this is the content of Scheme A. On the 
other hand, with the inclusion of Si mp in the calculation of stationary-phase points, the critical trajectories feature 
small deviations from retro-reflection, due to the difference in electron and hole wavevectors. (Note that Si mp is 
proportional to k + — /c_.) This, in turn, is Scheme B. For larger n (and thus higher resolution contributions to the 
DOS), Simp becomes larger, so that Scheme B should be used [Q. These higher-resolution contributions show up in 
the DOS as mesoscale oscillations due to the superconducting confinement of quasiparticles. 



A. Scheme A: Andreev approximation 



In the present section we shall evaluate Eq. (5.1) for the DOS within the stationary-phase approximation via 
asymptotic Scheme A. In this scheme, the stationary-phase points (i.e. the closed classical trajectories) are obtained 
by making stationary the phase SAnd alone. The factor expiS; mp is considered to be slowly- varying, and thus is 
evaluated at the critical points determined from SAnd alone. In all o ther factors, the difference between k+ and k- 



can be neglected. The reflection rule can be obtained from Eq. ( |4.6| ) by letting k + ,k^ — > kp. In this limit, velocity 
vectors are, upon reflection, exactly reversed. The classical trajectories obtained by this reflection rule are tracings 
of the chords of the billiard. This allows us to label every classical trajectory by two boundary points, along with 
the number of reflections (or tracings). Therefore, the closed classical trajectories of Scheme A with n tracings are 
degenerate (in the sense that their phases SAnd are identical) and they belong to a {2d — 2)-parameter family. 

Let us first explore the underlying physics of Scheme A. In conventional billiards, degeneracy of closed trajectories 
is usually related to the symmetry of the billiard (e.g. rotational symmetry of a circular billiard). However, in 
Andreev billiards this degeneracy is due to the underlying electron-hole symmetry at the Fermi surface, and is broken 
explicitly at nonzero energies by the term Si mp ]3^ ]. This (approximate) non-geometric symmetry is the reason for 
the (approximate) integrability of Andreev billiards, whatever the shape of the billiard. 

The method for evaluating the DOS is as follows. We fix two reflection points (so that the degeneracy is lifted), 
evaluate the rest of the (surface and volume) integrals in the stationary-phase approximation, and then evaluate the 
remaining two surface integrals. Without loss of generality, we choose the first and last reflection points as those to 
be fixed, and label them a and (3, respectively. We focus on the contribution to the DOS from the term having 2n + 2 
reflections: 

p 2n (E) = Im— — e 2mip Jdcr a dapd<j ai ■ ■ ■ da a2n dg T {ac,ac 1 )dg±{ai,at 2 ) ■ • • dg T {a 2ni P) J d d x g±(/3,x) <9g±(x, a), 

(5.7) 

where the +/— sign represents the contribution from the electron/hole sector. 

We first evaluate the x integral. The stationary-phase points of the x integration lie on the line joining a. to (3. 
Because the phase does not vary as one moves the point x along this line, it is natural to separate the x integration 
into longitudinal (i.e. parallel to a — (3) and transverse (i.e. perpendicular to a. — f3) components; see Fig. [?] for the 



nomenclature and geometry). We use the asymptotic expression (4.G) for the homogeneous Green functions at large 
argument, together with the asymptotic expression for dg, i.e., 

2 9<7±(/3,a)« {^-) 2 cosd «P exp(±ik±la0Tin{d-l)/4), (5.8) 

where 9 a p is the angle between the ray joining a to /3 and the normal direction at the point a. Then, we fix I, 
evaluate the t integral in the stationary-phase approximation, and evaluate the I integral: 



26 



FIG. 7. Geometry for x integration in Eq. (5.9) 
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(5.9) 



Next, we use Eqs. (5.8) and (5.9), together with Eq. (5.7), to obtain the following asymptotic expression for p 2n (E): 

%d—n—\ 



pf n (E) w Re J da a dap (cos8 a p cos6p a ) n J n (a,/3; kp 



( fcf 



V27r£ Q 



exp 



^— i2ri(y9 + in(k + — k-)t a p^ , (5.10a) 



/„(<*, /3; fc F ) = / da ai ■ ■ ■ da a2rl _ 2 

J sp 



<SAad 



(5.10b) 



and J indicates that the surface integrals over a.\ ■ ■ ■ a.2n-2 should be evaluated in the stationary phase approxima- 
tion. In order to do this evaluation, we expand the functions {tataA around (a2i,ct2i+i) = (ct,/3). 

We now focus on the evaluation of the stationary phase integrals in the definition of J n (o;, /3; /cf)- In the remainder 
of this section we shall work in three dimensions; after going through the calculation for this d = 3 case, the extension 
to higher dimensions is straightforward. The coordinate system for evaluating these stationary phase integrals is as 
follows (see Fig. ||). We use a pair of coordinate systems: one for the set {cn 2 i} (which are near a); the other for 
the set {a.2i+i} (which are near (3.) Each of these coordinate systems is constructed from the normal vector at the 
point in question, along with two unit vectors in the tangent plane at this point. Thus, the coordinate system doublet 
comprises two sets: one formed by the unit vectors {x Q ,y Q ,n Q }; the other formed by the unit vectors {x^, y^, fig}, 
where x a /^ and y a /p respectively lie on the tangent planes at the points ot/(3. The choices of orientation of the axes 
within the tangent planes are somewhat arbitrary; to ease the calculation we shall choose x Q to be parallel to xg 
(and thus parallel to the line of intersection of the two tangent planes), which fixes y a and yp. In this coordinate 
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FIG. 8. Geometry for fluctuations of the reflection points 



doublet, we shall label the local coordinates of point cti by (xi,yi). The three-dimensional position vectors are then 
determined as follows: (i) If i is even, use the coordinate system at point a; otherwise use the coordinate system 
at point p. (ii) Then the three-dimensional position vector in the chosen coordinate system can be expressed as 
8 t = (x i ,y l ,z Tlt (xi,y i )), where 

{a, for i even, , 
(3, for i odd, (5 - iij 

and z m (x, y) is the local equation of dV near the point rji. We first expand ( a2 ia 2 i+i to second order: 

e ~o rx x a i l^'-i ~ *^| 2 (n a)3 ■ S 2 j-i - n a0 ■ 8 2i ) 2 (k-\o»\ 
W-iQ 2i = «a/3 + n Q/3 • (d 2 i-i - d 2i ) H —„ —. , (5.12a) 

ZZ a p Zla/3 

f ~ V J-n IX X N , l^H!±l_^il!. (na/3 ■ $2i+l ~ n a/ g ■ S 2i ) 2 fKIOHl 
«Q2 I a 2 i + i = *Q/3 + n Q/3 ' (0 2i +l - 2i ) H — — , (5.12b) 

where n Q/ 3 = (/3 — a)/ 'tap. We then use this approximation to re- write ^And^ 

^ n— 1 

'S'Andrccv ~ -"J ( #2i-l ■ <$2i + (<5 2 i-l ■ np a )(S 2i ■ n a p) ) . (5.13) 

t-a/3 , V ' 

Note that z m (x,y) is a quadratic function of x and y and therefore, for the purposes of expanding SAnd to second 
order in {xj,?/i}, we can neglect the n I)i component of Thus, the curvature of dV does not feature in I n (a,f3; hp), 
which implies that the stability of these trajectories does not depend on the curvature p9| . Then, S^nd can De written 
as follows: 

^Andrecv ~ ~„ (^2i-1 2/2i-l ) ' D • ( 2 * ) , (5.14a) 

D _/l-n 2 -n x (n y cos0 + n 2 sin0) \ (5.14b) 
I — n x n y cos </> — n y (n y cos + n z sin 0) y ' 

where n x , and n z are, respectively, the x, y and 2 coordinates of n a| a) in the coordinate system at point a, and 
cos <j> is the angle between the y axes of the coordinate systems at the points a and (3. Then 

/2n-2 
JJ dxi d Vi exp ((ifc F /£ Q/3 )X t ■ V(n) ■ x) , (5.15a) 
»=i 
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Thus we arrive at the following expression for /„: 
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(5.15b) 



(5.16) 

where sgnl? denotes the signature of the matrix D (i.e. the number of positive eigenvalues minus the number of 
negative ones). By inserting this expression into Eq. (5.10a) we obtain the contribution to the oscillatory part of the 
DOS associated with the 2n-reflection term: 



pf n ( E ) ~ Re / da a da l 



hp COS 9 a p COS 9p a 

2£ a/ 3 



exp (— 2intp + m(fc+ — k)£ a p) . 



(5.17) 



By collecting together the contributions from all numbers of reflections (n = 1,2,...) we arrive at the following formula 
for the oscillatory part of the DOS in Scheme A: 



8p 1 (E) sa / da a dap — Re 



A similar calculation for the d-dimensional case leads to the result 
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(5.18) 



(5.19) 



Formula (5.19) for the oscillatory part of the DOS simplifies further in the 7-^0 limit. To see this, consider the 
following familiar identity: 
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(5.20) 



Applying this identity to Eq. (5. IE), we obtain 
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(5.21) 



The last term in Eq. (5.21) is a constant term, and it exactly cancels the leading-order Weyl (i.e. bulk) term. In 
order to see this, consider the coordinate transformation from (a, (3) to (b, n), where n is the direction of the chord 
and b is the position-vector specifying the intersection of the chord with the plane perpendicular to n (i.e. the impact 
parameter). The transformation of the surface elements are as follows (see Fig. |9|): 



da a dap cos 9 a p cos 0p a i 1 a p d = dndb. 



By using the (b, n) coordinate system, the last term in Eq. (5.21) can be cast into the following form: 
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(5.22) 



(5.23) 



A direct comparison with the leading-order Weyl term, given in Eq. (|5.2[ ), shows that the term in Eq. ( 5.23 ) identically 
cancels with the leading-order Weyl term. Thus, the DOS in Scheme A can be written as 



p(E)v £ 



k d ~ 2 f a 



m— — 00 



2(2tt) 



d-l 



<ln db () ( ~j^~£a0 — 2ip — 2wm 



(5.24) 



which is precisely the result found via the Andreev approximation, as stated in Sec. II C . 
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FIG. 9. Jacobian of the transformation from boundary points to scattering parameters 



B. Scheme B: Mesoscale oscillations beyond the resolution of the Andreev approximation 

The main motivation of Scheme B is to capture the mesoscale oscillations in the DOS that are caused by confinement 
of quasiparticles by the superconducting surround. The reason that Scheme A (and thus the Andreev approximation) 
is not capable of capturing such oscillations is that in Scheme A the transverse degrees of freedom are not quantized 
(i.e. quasiparticle motion on the chords is quantal, but there are chords arbitrarily close to one another, indicating 
that the transverse degrees of freedom are treated classically). On the other hand, in Scheme B we take into account: 
(i) the imperfectness in retro-reflection (arising from the previously-neglected difference between the wave vectors of 
incident and reflected electrons and holes); and (ii) the imperfectness in charge-interconversion, the amplitude for 
which is C(A//i). A priori, we know that mesoscale oscillations in the DOS must originate from one or more of the 
processes not yet taken into account, these being (i) and (ii), above, as well as the subleading quantal corrections 
(that are ignored in both schemes). In fact, as we shall see, it is the imperfectness in retro-reflection that is primarily 
responsible for the mesoscale oscillations in the DOS. (Note that the imperfectness in retro-reflection occurs transverse 
to the incoming direction.) The imperfectness in charge-interconversion does modify these DOS oscillations. However, 
by itself (such as in a model with perfect retro-reflection but imperfect charge-interconversion) it is not capable of 
producing them. In order to clarify this issue, and thus to assess the significance of these processes, we shall define a 
useful intermediate model: the Perfectly Charge-Interconverting Model (PCIM). This PCIM has the feature of being 
fully quantum- mechanical; however, in it, any single reflection from the boundary is certain to have the effect of 
converting electrons to holes (and vice versa). As it contains all quantal effects, the comparison of its predictions with 
those from the semiclassical approach enable us to assess the importance of quantal effects beyond the semiclassical 
limit. Moreover, from this comparison it is possible to draw conclusions regarding whether it is imperfectness in 
retro-reflection that is capable of capturing the mesoscale oscillations in the DOS. 



1. Perfectly Charge-Interconverting Model 

We are now at a position t o def ine a Perfectly Charge-Interconverting Model (PCIM). We start with the expansion 
for G in terms of G R in Eq. ( 4.50| ). However, we shall replace G R by its leading-order form, i.e., — ie _l ^o"iG N . Then 



the resulting model is defined as an integral equation for G, residing inside the billiard: 

G(x,x') = G N (x,x') - 2ie- lv [ da a <9G N (x, a) a x G{a, x'). (5.25) 

JdV 

The off-diagonal matrix <j\ ensures that, upon each reflection from the boundary, electrons are totally converted into 
holes (and vice versa). Moreover, this model is fully quantum-mechanical, in the sense that it retains wave-propagation 
effects (as implied by the presence of surface integrals.) Let us now focus on the DOS of the PCIM, treated at the 
semiclassical level. 
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Due to the imperfectness in retro-reflection in Scheme B, the corresponding classical dynamics is no longer a priori 
integrable; on the contrary, it is weakly chaotic for generic shapes |2j. However, the closed periodic orbits do fall 
into two quite distinct classes: one consists of multiple tracings of each stationary chord (i.e. chord of stationary 
length, which we refer to as a E); the other of much longer trajectories that "creep" around the billiard boundary (see 
Fig. ^J). Correspondingly, the DOS is the sum of: (i) an average term, which depends on the volume of the billiard 
(i.e. the leading Weyl term); (ii) a finer-resolution term, having a universal lineshape that depends solely on the length 
and endpoint-curvatures of the Es; and (iii) highest resolution terms, which depend on the classical dynamics of the 
billiard in question. 



2. Stationary chords 

Here, we focus on the contribution to the DOS coming from stationary chords. As we shall derive below, it is 
possible to give a closed-form expression for the contribution of a E to the DOS for billiards of generic shape. It is, 
however, necessary to distinguish between isolated Es and degenerate Es. First, let us consider an isolated E in a 2D 
billiard, with endpoint curvatures R\ and R 2 . Then the E contribution to the DOS from 2n reflections is given by 

^ ± = ^-Rey-J (— j exp^-M^-^)/. (5.26) 
Here, is the Gaussian integral resulting from the expansion of the action to second order: 



f n { k k 

In = J \\_dxi exp^ i ( -j-xn-x x 2l - -j^x 2i x 2i+1 



where X2n+i = x\. In order to evaluate this integral, we define a matrix M such that: 

f k + k- k + — fc_ / 1 1 \ , / 1 1 \ , \ 
Xi M t j Xj = ^ » ( J^ x K-i x 2i - j^x 2i x 2i+ i + I — - — 1 x 2l _ x + I — - — ) x 2i ) . (5.28) 



i,j i=l 



(Note that, by definition, an isolated stationary chord is one for which none of the eigenvalues of M is zero; whenever 
a zero eigenvalue occurs, the stationary chord is said to be degenerate.) Then, for an isolated E , /„ can be expressed 
in terms of the determinant and signature of M: 

I n = ; = exp (i-senAl) . (5.29) 
v/|dctM| V 4 S J { ' 

The eigenvalues of M can be obtained from the eigenvalue equation . My Xj = Xxi . The symmetry of M under the 
transformation x% — > x i+2 restricts the form of the eigenvectors to be 

_ A i2m7r(|i) tj i2mir ( ) 

x 2 j — A m e \2-ni ^ x 2j+ i — ±S m e \ 2 ™ >, 
where m — 1, 2, ■ • ■ , n. Then the eigenvalue equation reduces to 



fc_i_— fe_ I \ 1 \ fe+ e _i w m fc_ 



L e -™? _|_ ^± e "r? fc+^fc- ( 1 1 \ / \ B m I m \ Br, 



2 \£s Ri 

The determinant of M is readily obtained as 

k- I 1 1 \ k+ .-j,™ k 



(5.30) 




{is R2 J 



ie - " 1 ^ - ^e" 



t\ — isRi — £t,R2 \ ( k + — k- \ 2 .k + k^ . 2 ( _m 



1 k±„in— k + —k- (11 

' A— e n i b_ e n n 



, — 4 — ^ — sin ( 7r 

-R1-R2 \ (t, iy V n 



(5.31) 



For n <g; yJk+kZ / '{k + — fc_), this expression can be simplified to read 
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detM 



4 - feRi - £ s i? 2 



R1R2 

( _ 1)n _l - fegl - fei?2 



f 2 



n -■(-=) 



ly 



k+k. 
1 s 



m— 1 
n— 1 



(5.32) 



The signature of M is obtained from the eigenvalue problem (5.30) via the inves tigat ion of the signs of the eigenvalues 
A^ [or, equivalently, from the trace and the determinant of the matrix in Eq. ( 5.3 Cl| ) ] . If the determinant is negative 
then A+ and A~ are of differing signs and, hence, the associated pair cancel in the signature of M. If the determinant 
is positive then there is a pair of negative or of positive eigenvalues, the sign of which is determined by the trace. 
These considerations allow us to express sgn M as follows: 



sgnM = J2 ( 1 



sgn 



ly - lyRi - lyRi 
R1R2 



4^ — sin 



I 2 

11 y 



■ 2 ( m\ 
m V n ) 



Again there is simplification for n <C \J k + k-/(k + — fc_), because in this regime only m — n term contributes to the 
sum above and, hence, the expression for sgnM becomes 



sgn M w 1 + sgn 



R\ — i?2 



R1R2 



sgn 



*E Rl 



1 

~R2 



In particular, when R\, R2 > we have 



sgnM w sgn (R 1 + R 2 - ly,) - 1. 



(5.34) 



(5.35) 



Putting all the pieces together, the expression for the contribution to the DOS originating from isolated stationary 
chords is as follows: 



Py.± 



27rfc ± Re UJ \iy 



n/2 



k. 



n/2 



1 



yj\ detM 



exp 



(in(k + — k-)£y — i2nip + i— sigM^j 



(5.36a) 




k + k^£sRiR 2 



4Tr 2 (k + -k-) 2 \£y-R 1 -R 2 



exp (in(k+ - k_)ly - i2ntp + i^(sgn (R 1 + R 2 - ly) - lfj, (5.36b) 



where the second line is valid for n <C ^Jk+kZ j '(k+ —k~). We can now sum this expression to all orders in n to obtain 
the general form for the contribution to the DOS originating from an isolated stationary chord: 



Re- 



(fc 4 



-) 2 lyRiR2 



4Tr 2 k + k-(k + - k_) 2 \ly -Ri-R 2 \ 



if (8gn(Ri+-R 2 -fe)-l) 



ln(l 



i(fe_l_ — —i2rnp 



(5.37) 



The DOS expression in Eq. ( [5.37 ) is valid when its prefactor is free of singularities. Such singularities would indicate 
that the fluctuation determinant has a vanishing eigenvalue. This would mean that the stationary chord is no longer 
isolated; instead there is a direction, corresponding to the vanishing eigenvalue, along which fluctuations do not change 
the phase. 

Let us now focus on these prefactor singularities. The first kind occurs when k + = fc_ or, equivalently, when 
E = 0. The reason that this singularity arises is that when k+ = fc_ retro-reflection is perfect and, thus, all chords 
become stationary. As, in Scheme B, we are implicitly assuming that fc + and fc_ are distinct, this kind of singularity 
is an artifact of the approximation scheme, and thus is unphysical. The second kind of singularity occurs when 
ly = R% + R 2 . This is a geometrical singularity, in the sense that whenever the shape of the billiard is such that 
chords of stationary lengths are not isolated, such a singularity arises. For example, if the billiard is circular, the 
chord of stationary length is the diameter of the circle, and its length does not change when it is rotated. In this case 
R\= R 2 = ly/2. Another example is shown in Fig. [10| see the caption for details. Such singularities can be handled 
in Scheme B as follows: we fix one (or more, if necessary) reflection points, so that the £ is no longer degenerate; we 
then integrate the remaining surface integrals in the stationary phase approximations; and, after that, we evaluate 
the remaining surface integral. In this way, we obtain the general form for the stationary-chord contribution to the 
DOS: 



py(E)^RcJ2Zy c 4A7r / 4 L 



id-i- 



-(1 



fe_ — 2itp 



(5.38) 
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FIG. 10. Example of a billiard that is not circular but has degenerate stationary chords. Opposite segments of the boundary 
are partially coincident with a pair of concentric circles, as shown; when this happens, Is — Ri + 7?2- 



Here, Li„(z) = J2JLi zJ 7j™ 1S the polylogarithm function, w is the dimensionality of the degeneracy of the E (e.g. w = 1 
for a circle), is a slowly varying real function of energy, which determines the size of the DOS oscillations, and 
A is a measure of the stability of the Es, which determines whether the "tail" goes to higher or lower energies. For 
example, for a degenerate E and d = 2 we have 

/ (fc+ + k.)HvR 
S S Y 87r 3 fc + fc_(fc+-fc_)|£ s -i?|' 

A = sgn(i? — £y), 
where Vy, is the volume of the degenerate E and R = (R% + R-i)j2. 



3. Creeping orbits 

In this section, we shall obtain the finer oscillatory structure that lies beyond the stationary-chord contribution. 
The method for obtaining this structure consists of: (i) finding the classical periodic orbits; (ii) evaluating the surface 
integrals in the stationary phase approximation (i.e. expanding the action to second order around each classical 
periodic orbit, thus reducing the surface integrals to Gaussian integrals, and then evaluating the resulting Gaussian 
integrals); and (iii) summing over all periodic-orbit contributions. For the purposes of illustrating this method, we 
now focus on the case of a circular Andreev billiard and obtain an expression for this finer oscillatory structure in 
DOS. In addition to its illustrative purposes, having obtained this expression for the finer DOS oscillations for circular 
Andreev billiard will become useful when we discuss an additional approximation, valid for shapes with slowly varying 
curvatures, as we shall do later in this Paper. 

(i) Classical Dynamics: The classical dynamics of Andreev billiards is defined through the Scheme B reflection rule: 
k + sin6> + = fc_ sin6>_. Here, #+/_ is the angle of incidence (or reflection) for the electron/hole. In the case of circular 
Andreev billiards, the classical periodic orbits can be specified (up to rigid rotations of the full trajectory) by the 
number of reflections (2n) and by the number of times the orbit "winds" around the billiard (j). The fact that 
classical periodic orbits of circular Andreev billiards can be specified by (n,j) is a consequence of the integrability of 
the classical dynamics. 

Our goal is to express dynamical quantities, such as the action S c \(n, m, fj,, E, R), angles of incidence and reflection, 
and the length of propagation between two successive reflections (which we shall denote by £±) in terms of the 
parameters (n,m, fi, E, R). Here, R is the radius of the circular billiard. First note that, for a given orbit, the angle 
at which an electron is incident or reflected is the same at all reflections; the same fact holds for holes (see Fig. pT[ ). 
Next, note that each reflection rotates the direction of the momentum by 69 = 6- — 9 + . Thus, in order for the orbit 
to close after In reflections the momentum direction must be rotated by 2-7rj. In other words, SO = n j jn. However, 
6± , must also satisfy the classical reflection rule k+ sin 8 + = fc_ sin . By using these conditions we find (after a little 
algebra) that for a closed orbit specified by n and j we have 

cosg ± ^± fc± ~^ A C » s( " j/n) , e ± = 2Rcos6 ± , (5.39) 
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FIG. 11. A closed periodic orbit with imperfect retro-reflection; the degree of imperfectness is exaggerated for the purpose 
of illustration. 



/ \l/2 

where Afc = I /j, — 2^/ /i 2 — E 2 cos(nj/n) J . Observe that, owing to the fact that cos6* + > cos0_ > 0, the possible 
values of n and j are restricted: 

cos (l[f)>^. (5.40) 
\ n k. 



Next we evaluate the action: 



S c \(n^ m, /it, E, R) = nk^i^ — nk—£- 

= 2n(k + Rcos6 + — fc_i?cos0_) 

1/2 



= 2nR Ut - 2-v/m 2 - E 2 cos(nj/n)^ 
= 2nAkR. (5.41) 



(ii) Evaluating surface integrals: Having obtained the critical points of the surface integrals in Eq. (5.1 



proceed to expand the action around these points. In order to do this, we expand each g±(a,/3) in Eq. (5.1) around 
the critical values of its arguments (i.e. (a c ,/3 c )). This can be accomplished by writing 



a = at c + d a , (3 = (3 C 



(5.42) 



and expressing 8 a /p in the coordinate system at at./ (3 with coordinate axes specified by the normal (i.e. n a ^) and 
tangent (i.e. t a //3) directions of dV at a/j3. In this coordinate system, the coordinates of the boundary point 8 a /p 
can, up to quadratic order, be parametrized as 



*//3 



2R r 



IP 



(5.43) 



where R a /p are the radii of curvature at a./ (3. The surface elements can also be expressed using the parameter s a /p: 

da a//3 = ds a//3 + 0(l/R). (5.44) 

Then dg± can be approximated as 
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FIG. 12. A closed periodic orbit with imperfect retro-reflection; the degree of imperfectness is exaggerated for the purpose 
of illustration, but less so than in Fig. [lT| 



2%fc(s ct ,5 / 3) ~ Ty ^ exp (±ik±e ac p o =F iTr/4±ik±$(s a ,sp)^, (5.45) 

where the phase function $ is given by 

sj COS s 2 cos 9 a 1 
$(s a ,s,a) = sin 0£ - s a sin 9 a ^— h — (s/3 cos 6/3 + s a cos a ) . (5.46) 

ZHp ZU a Zi a< . p c 

The terms linear in s a /p cancel with the linear terms of the next and previous Green function, due to the stationarity 
feature of the critical point. For the circle we have a = dp = 6±, R a = Rp = R, and i ac /3 c — 2Rcos9±. This allows 
us to write 

cos 9 A- 

$±{s a ,sp) = —{s a -sp) 2 , (5.47) 

and thus 2dg±{s a , sp) — 2dg±(s a — sp). We are now in a position to evaluate the DOS oscillations due to a periodic 
orbit with 2n reflections that winds j times around the circular billiard: 

p± = 2 ll k± Rc exp y Scl ( n ii> ^' E > R } _ i2ni Pj In , (5.48) 

- 2n 

I n = n\ds i 2dg + {s 1 )2dg-{s l -s 2 )---2dg + {s 2n - l ). (5.49) 

J i=l 

Although it is possible to evaluate this Gaussian integral by the usual formula that relate s it to the determinant and 
signature of the quadratic form in the phase, the (rotationally invariant) form of <f> in Eq. Q5.47 ) allows the evaluation 



of the integral above in an easier way, viz., via Fourier decomposition. By defining the Fourier transform 

/oo 
dse ips 2dg ± (s) = iexp (±ip 2 R/k± cos6» ± ) (5.50) 
-OO 



we are able to write I n j as 



/,,; - / ^-(-2dg + (p)2dg-(p)y 
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/ 7T- (-iT cxp (inp 2 R ({k+ cos 0+y 1 - (k- cos 6L)- 1 
J-oo 27r 

(-0" 



k + cos + fc_ cos 



AnnR(k + cos + — A;_ cos #_) 



exp —in /A. 



(5.51) 



Then the contribution to the DOS oscillations originating from a creeping orbit with 2n reflections and winding 
number j can be written as 



«j 



Re 



' cos 9 a 



R 3 k + cos 9 + k- cos 6>_ 
nn(k+ COS0+ — /j_ cos#_ 



-(-i)"exp (^^(n, j, ^i, i?) - 2inip - in/Aj 



= Re 



'i? 3 (fc+ — fc_ cos(nj/n)) (k + cos(nj/n) — fc_) 



nn (A&) 5 



E cos(nj /n) exp (^2in Ak R — 2intp — in / A 



(5.52) 



(wi) Summing over all creeping orbits: In order to ob tain the creeping-orbit contribution to the DOS oscill ations , 
we must sum the contributions specified by Eq. (5.52) for all n and j, obeying the restriction given in Eq. (5.40). 
Although this summation may appear to be problematic, we remind the reader that E has a nonzero imaginary part 
r, viz., the smoothing width, which suppresses exponentially (in n) the contributions coming from high n values. 
Thus, for a given T it is possible to perform the sum up to a value of n such that any contribution from higher n 
would be invisible on the scale set by the truncated sum. 
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FIG. 13. Density of states oscillations for a circular Andreev billiard: 

r/n = l.i x icr 4 . 



k F R = 150; A/u = 0.08; smoothing width 



In Fig. |l^ we plot three versions of the oscillatory part of the DOS: the Scheme A result (dashed line), which includes 
all chords but is dominated by the stationary ones; the Scheme B result (dotted line), which includes creeping orbits 
and stationary chords; and the exact PCIM result (full line), obtained by the numerical solution of the PCIM. Observe 
that, although the local average behavior of the exact DOS is essentially that captured by Scheme A, in order to 
capture the mesoscale oscillations beyond this average behavior one must use Scheme B. 

Poisson summation and the semiclassical quantization condition: In the case of a circular Andreev billiard the 
summation over n and j can be performed approximately to all orders. The procedure for doing this involves using 
the Poisson summation formula and evaluating the integrals in the stationary phase approximation. By using this 
procedure it is possible to obtain the energies at which the DOS has simple poles, viz., it is possible to obtain a 
semiclassical quantization condition. Our starting point is the expression Eq. ( |5.52 ) for p n 'K In terms of this, the 
DOS oscillations can be written in the form 



Rc 



a(nj/n) 



■ exp 



(i2nAk(nj /n)R — 2inip — in j '4 



(5.53) 



where the amplitude a is defined via 



a(x) 



/ i? 3 (/c + — k^ cos(a;)) cos(a;) — fc_) 
= V nAk(x) 5 k+k- 



Ecos(x). 
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TABLE I. Comparison of the eigenenergies of the PCIM for a circular Andreev billiard computed semiclassically and com- 
puted exactly (numerically) for kpR = 150; A/p = 0.08. 



m 


n 






AE/p 





4 


0.066694 


0.066695 


0.000001 


1 


4 


0.066698 


0.066696 


0.000002 


2 


4 


0.066699 


0.066700 


-0.000001 


3 


4 


0.066708 


0.066707 


0.000001 


4 


4 


0.066716 


0.066716 


0.000000 


70 


4 


0.073978 


0.073967 


0.000011 


71 


4 


0.074204 


0.074206 


-0.000002 


72 


4 


0.074444 


0.074450 


-0.000006 


73 


4 


0.074711 


074700 


00001 1 


74 


4 


0.074951 


0.074954 


-0.000003 


110 


3 


0.067310 


0.067285 


0.000025 


111 


3 


0.067842 


0.067880 


-0.000038 


112 


3 


0.068466 


0.068495 


-0.000029 


113 


3 


0.069169 


0.069132 


0.000037 


114 


3 


0.069840 


0.069791 


0.000049 



a PCIM 
b Semiclassical 



By using the Poisson summation formula we first cast the expression for the DOS in the form 

dj a[7TJ l_ n) exp (2mAfc(7rj/ n)R — 2intf — 2inmj — i7r/4) . 

Here, m has an interpretation as the angular momentum quantum number. We then evaluate the j integral in the 
stationary phase approximation. The stationary phase points satisfy 



= f(m,R,p,E), (5.54) 

where A = m 2 /(k^R) 2 , and only real and positive values of cos(7rj c /n) are allowed. The important point to observe 
here is that cos(irj c /n) is independent of n, which implies that Ak is independent of n, too. This allows us to write 
the expression for p as 

p w {^k F R{f{m) - X)y 1/2 a (cos' 1 f(m)) Re ^ e m ( 2Afe(m) ^ 2 * , - 2mcos_1 f ^). (5.55) 

tci n 

We are now in a position to perform the sum over n: 

i(2Afe(m) fl-2ip-2mcos _1 f{m)) 

p^Y, (irk F R(f(m) - A))" 1 2 a (cos" 1 f(m)) Re - - e<(2Afc(tn) fl _ 2y _ 2mcog -i /(m)) ■ (5-56) 

m 

Thus, the semiclassical approximation to the eigenenergies for a given angular momentum quantum number m (i.e. the 
DOS peaks) are given implicitly as the roots of 

expi(2Afc(m)i?-2^-2mcos" 1 /(m)) = 1. (5.57) 

In Table | we compare the eigenenergies of the PCIM calculated via this semiclassical scheme and exactly. From this 
Table we see that, as expected, the semiclassical results agree with the exact results upto contributions of relative 
order 1/kpR. 
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4- Incorporating ordinary reflection 



Thus far in our semiclassical treatment we have ignored all amplitudes involving ordinary reflection. For non- 
grazing incidence [i.e. for 9—(jr/2) ~ 1] the amplitude for ordinary reflection is very small (in fact, of order A//2 cos 2 9). 
However, orbits that contribute dominantly to the oscillatory structure of the DOS obey \0—(w/2)\ <C 1, and the refore 
ordinary reflection amplitudes are not negligible and must be incorporated. This can be done by returning to Eq. ( 4.50| ) 



and re-evaluating the trace formula using the full expression for G R (i.e. not just the leading, off-diagonal, term). 
Then, as a result of treating both charge-interconverting and charge-preserving reflections, the classical limit changes 
drastically: the initial conditions of specifying position and momentum (and charge, in the case of Andreev billiards) 
no longer determines the full orbit; on the contrary, each reflection splits the incoming ray into two rays: one (albeit 
imperfectly) retro-reflecting; the other, ordinarily reflecting. Thus, the classical dynamics is no longer deterministic 
(i.e. specifying the position and momentum no longer determines the full orbit) and, instead, the initial conditions 
specify a superposition of orbits. (A similar situation emerges in the context of Schrodinger billiards when the billiard 
has sharp jumps in the single-particle potential |jt6|.) Generically, the number of closed orbits increases exponentially 
as a function of number of reflections (as opposed to linearly, as is the case when ordinary reflection is neglected) , but 
the amplitudes for these orbits are suppressed exponentially (owing to the fact that the amplitudes for ordinary and 
Andreev reflections are smaller than unity), allowing the number of reflections to take care of any divergences arising 
from this ray-splitting feature. 

The exponential increase in the number of closed orbits, as a function of the number of reflections, makes evaluation 
of the periodic-orbit sum difficult. For sufficiently smooth shapes there is a way to circumvent this difficulty, which 
involves resorting to a different approximation scheme, as we shall shortly show. The main motivation for this 
approximation scheme is as follows. 

1. The closer a primitive creeping orbit is to the boundary (i.e. the closer \6— {tt/2)\ is to zero), the shorter it its total 
length. As the periodic-orbit contributions are suppressed exponentially with their lengths, the creeping orbits that 
are closer to the boundary contribute more strongly to the DOS oscillations. 

2. The closer a creeping orbit is to the boundary, the bigger the ordinary reflection amplitude (~ A/^cos 2 6.) Thus, 
the orbits that involve ordinary reflection contribute more strongly if they are close to the boundary. 

3. For the orbits close to the boundary (which, by the previous considerations, dominate), consecutive reflections take 
place very near to each other, and thus "see" only the local curvature of the boundary. 



These considerations motivate us to perform an "adiabatic" approximation to the expansion in Eq. (4.50), in which 
we assume that the curvature of the boundary varies slowly, relative to the rate at which creeping orbits sample the 
boundary. 

Our starting point is the integral equation ( 4.37| ) that generates the MRE: 



G u (x, x') = G N (x, x') + f da a d a G N (x, a) • /x u (a, x'), 

JdV 

H il sa 2G R + 28G K n u , (5.58) 



where G R is given in Eq. (4.41a). These equations can be cast to the following form: 

G u = G N + 2<9G N (f - 2<9G R ) _1 G R , (5.59) 

which shows that the approximate poles of G 11 are given by the poles of K = (l — 29G R ) . Thus, in order to obtain 
the energy eigenvalues it is sufficient to consider the following integral equation defined on the surface dV: 



K(a,/3) =IS(a,P) + 2 [ da 1 <9G R (a,7) K( 7 ,/3) 

JdV 



This equation will have a regular solution (i.e. one without poles) if and only if none of the eigenvalues of the operator 
2<9G R is equal to unity; conversely, poles of K occurs at energies for which at least one of the eigenvalues of 2<9G R 
is equal to unity. Consequentially, we now focus on the following eigenvalue problem defined on the boundary of the 
billiard: 

2 I do- dG R (a,/3)u(/3) = Au(a). (5.60) 

JdV 

We shall work in the coordinate system in which the boundary is parametrized by its arc length s. (Recall that we 
are considering 2D billiards.) Thus, the equation for the boundary is given by the vector function a(s), 
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t(s) = da/ds 

is the tangent vector, 

n(s) ee -R(s)d 2 a/ds 2 

is the unit normal vector, and 

R(s) ee d 2 a/ds 2 1 

is the curvature at the point ot(s). Furthermore, we shall transform to the coordinates 



and 



S = 



and define 



G (t, S) ee G R (s,s') 



(5.61) 



(5.62) 



(5.63) 



(5.64) 



(5.65) 



The virtue of this coordinate system is that for a constant-curvature boundary (viz. a circle) G R (s, s') = G R (s — s'). 
Thus, if the curvature is slowly varying then G R (t, S) is a slowly varying function of S. In this coordinate system, 
then, the eigenvalue equation Eq. Q5.60p becomes 



dtdG K (t,s-t/2)u(s-t) = Au(s), 
u(0) = u(L), 



(5.66a) 
(5.66b) 



where L is the length of the boundary. From periodic-orbit theory we already know that the dominant mesoscale 
oscillations in the DOS arise from the part of phase space in which the component of the excitation momentum lying 
tangent to the boundary is O (k-p). Thus, in order to capture the dominant mesoscale oscillations, it is sufficient to 
solve the eigenvalue equation ( 5.66a ) for the sector of eigenfunctions varying on the length scale 0(kp). In this sector 
of rapidly-varying eigenfunctions, only the small-i behavior of the kernel dCx(t, S) is relevant. In the following, we 
shall obtain an expression for dG (t, S) valid for small t. 

Our starting point is the general expression for G in polar coordinates, 



G(s, S ')=E R ' 



| J m (k + r < )H m {k + r > ) 





j-J m (fc_r<)iJ m (fc_r>) 



AraQ 



(5.67) 



where R(m/i?(s)) is defined in Eq. (4.41b). For s near s' (i.e. for small t), one can choose the "polar" coordinate 
system, in whic h a circle of radius R((s + s')/2J coincides locally with the surface at point (s + s')/2. Then the 
expansion (5.67) can be written as 



G(t,S)^J2 K - 



\J m {k+R{S))H m {k+R{S)) 






{ J m (k^R{S))H m (k^R{S)) 



imt/R(S) 



(5.68) 



The normal derivative dG is given by 

-^(j m (k + R)H m {k + R)) 



0G(t,S)«£ R -7 



R=R(S) 



■£t(J m (k-R)H m (k-R)) 



Jmt/R(S) 



(5.69) 



R=R(S) 



Having obtained the approximate form of dG, valid for small t, we shall now seek the eigenfunctions. We shall assume 
that the rapidly- varying eigenfunctions have the following form: 



u(s) = u(s) exp ims/R(s), 
viz., a slowly- varying envelope and a rapidly varying part. Then the eigenvalue equation becomes 



(5.70) 
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Au( s )e ,m(s)/fl(s) =2 / dtdG R (t,s~t/2)u(s-t)e vm{s - t)/R{s - t \ 



imt/ R{s 



u(s), 



2ttR(s)-L 



Au(s) pb 2 / dtdG H {t,s)e 
Jo 

A(m, s) u(s) = 2 / <ft dG R (i, s) e - imt /-R( s ) u(s) 
./o 

/ .2ttH(s) , 

A(m, s) u(s) = V 2dG R (m', s) / di e"^™"" 1 )t/i?(s) - 25G R (m', s) / 

Jo Jo 

\(m,s)u(s)=J22dG R {m\ S )2TrR{s)6 mm ,-2dG R {m\s) / dt e -i<v>-™')*W>) q(s), 

pb 51 2dG R (m', s) (2TrR(s)5 m , m , - (2nR(s) - L)) u(s), 

A(m, s) u(s) « 2<9G R (m, s) 2tt.R(s) u(s) + O {\dR/ds\ L) . 
By ignoring the term O (\dR/ds\ L) we have reduced the eigenvalue equation to 



(5.71a) 
(5.71b) 
(5.71c) 

dte -i{m- m ')t/R{s) Q ^ ( 5 .71d) 



(5.71e) 
(5.71f) 
(5-71g) 
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u(s) = A(m, s) u(s), 



(5.72) 



The important point to note about Eq. ( 5.72| ) is that it is simply a 2 x 2 matrix eigenvalue equation and, thus, its 
solution is straightforward. After obtaining the eigenvalues A(m, s), the DOS is given by the following formula: 



p(E) 



L ds 



o *vR(s 



d\(m, s; E) 



dE 



(5.73) 
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FIG. 14. Density of states oscillations for a circular Andreev Billiard: kpR — 150; A/fi = 0.08; smoothing width 
T/n = 1.1 x 10" 4 . 



In Fig. [li] the DOS arising from the solution of the eigenvalue equation Eq. ( 5.72| ) is compared to the the one 
arising from the exact solution of the full BDG cigcnproblem for the case of a circular Andreev billiard, which we 
have computed numerically. 

In this section we have shown that mesoscale oscillations in DOS essentially arise from imperfectness in retro- 
reflection. However, in order to correctly account for these oscillations it is necessary to account for the effects of 
imperfect charge-interconversion as well. The latter effects can be incorporated via an extension of the trace formula, 
in which a generic closed periodic orbit has both charge-interconverting and charge-preserving reflections. However, 
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for billiard shapes having slowly-varying curvatures, it is possible to obtain an adiabatic approximation to the DOS 
that bypasses the periodic-orbit summation (which has a number of terms that increases exponentially with the 
number of reflections) and, hence, reduces the task to the solution of 2 x 2 matrix eigenvalue equation. The strength 
of this method is that it relates the DOS to the closed form of R(m/i?(s)) and, thus, it is readily extendible to cases 
in which the N-S boundary is not clean (i.e. reflection amplitudes are modified). As seen above, in order to calculate 
the DOS one simply needs the reflection amplitudes at the billiard boundary. 



VI. CONCLUDING REMARKS; PERSPECTIVES 



In the present Paper we have explored semiclassical approaches to the oscillatory part of the density of states of 
Andreev billiards. We have done this by deriving two semiclassical trace formulas, each corresponding to one of two 
limiting schemes of the physical para meters specifying the billiard. The first of these trace formulas (viz. the Scheme A 



trace formula, discussed in Sec. V A) is essentially equivalent to the conventional quasiclassical approximation scheme 
first introduced by Andreev. The physical ingredients of this scheme are perfect charge-interconversion and perfect 
retro-reflection. It captures th e co arsest oscillations in the DOS. The second trace formula (viz. the Scheme B 



trace formula, discussed in Sec. VB ) not only captures the coarsest oscillation, but also goes beyond this resolution 
by capturing mesoscale oscillations. At the semiclassical level, mesoscale oscillations arise from orbits featuring 
imperfect retro-reflection. Although such oscillations are sensitive to charge-preserving reflection amplitudes (in 
addition to charge-interconverting reflection amplitudes) from the N-S interface, they are present even if there is no 
charge-preserving reflection. 

The methods developed in the present Paper readily apply to settings such as the superconducting proximity effect 
or the Josephson effects. Cases in which the phase of the pair-potential is relevant can be addressed by the appropriate 
mod ificati on of the renormalized Green function for complex A. In particular, by using t he m ethods described in 



Sec. IVD, it is possible to show that the leading-order behavior of G R is modified [cf. Eq. (4.49)], becoming 



/ o -g-(a x')\ 

G>,x')«exp(-^ + ^) ' , (6.1) 

Vs?(a,xO / 

where the phase <p = cos _1 (_E/|A|) and the phase 4> is the phase of the pair-potential. Thus, as one might have 
anticipated, the phase of the pair-potential simply adds to the phase acquired by the reflection. By using this modified 
form for G R one can, e.g., account for the zero-energy states observed in 7r-junctions (i.e. Josephson junctions through 
which the pair-potential undergoes a single sign-change). 

The Multiple Scattering Expansion that we have developed applies to systems consisting of piecewise homogeneous 
N or S regions. It is possible to generalize this expansion to handle smoothly- varying A(x). This can be achieved via 
an "energy-slicing" construction (rather than the time-slicing kind familiar from, say, the derivation of the Feynman 
path integral). In this way, one arrives at a functional version of the Multiple Scattering Expansion. To get a feeling 
why, let us divide the pair-potential range (0, Ao) [where Ao = maxA(x)] into a large number N of equally-spaced 
subintervals. Let us also break the x space into subregions, in each of which the pair-potential has values lying in 
solely one of the N energy subintervals. For N large, A(x) can be taken to be constant in each subregion. Then 
we can apply the Multiple Scattering Expansion formalism for this intermediate system to obtain its Bogoliubov- 
de Gennes Green function. Then, by taking the limit N — > oo, one recovers the full Bogoliubov-de Gcnnes Green 
function associated with A(x). Other spatial inhomogeneities, such as in the single-particle potential, can also be 
handled this way. Semiclassical approximations, as discussed here, can be applied to this expansion. For systems 
that have distinct regions of N and S, but in which these regions are modified from the ideal piecewise-homogeneous 
state because they feature a smooth variation of pair-potential, can still be regarded as Andreev billiards. In fact, 
as the amplitude for charge-preserving reflection is expected to diminish for more slowly varying pair-potentials, we 
expect the Perfectly Charge-interconverting Model (PCIM) to be more accurate for such systems. On the other hand, 
billiards constructed by fabricating a normal region inside a superconductor should feature more charge-preserving 
reflection, owing to interface effects. Such cases can be modelled by suitably modifying the reflection matrix R in the 
definition of the renormalized Green function G R . 

Apart from problems related to superconductivity, the methods presented here are also applicable to more con- 
ventional topics in quantum chaos. One such application is to the so-called Ray-splitting billiards. These billiards 
consist of regions of piecewise homogeneous (single-particle) potential, and the (sharp) boundaries between these 
regions serve as ray-splitting boundaries. By changing the homogeneous N and S Green functions to the homogeneous 
Helmholtz Green functions appropriate for a given constant potential, the formulation in the present Paper is readily 
extended wM- Another potential application is to multidimensional tunneling, studied in path- integral language in 
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Rcf. Ml ■ The Green- function language adopted in the present Paper is especially well-suited to the study of tunneling, 
owing to the energy (rather than time) dependence of the Green function. 

One experimentally relevant application of this work is to antidot Andreev billiards. Such billiards (in particular 
their magnetotransport properties) were recently studied experimentally by Eroms et al. |B2| . In this realization of 
Andreev billiards, it is the S region that is embedded in the N region, rather than the converse. (An experimental 
virtue of this geometry is that it is well suited to the study of the effects of weak magnetic fields, as the magnetic flux 
need not be quantized.) The methods presented in th present Paper are readily applicable to this antidot geometry. 
In particular, in App. |C], we describe the new features that emerge when the N region is nonconvex, as it is in antidot 
billiards. Moreover, the incorporation of magnetic fields-at least in the weak field case, so that it is simply excluded 
from the S region-is handled by modifying the Green functions (G N and G s ) to include the magnetic field. As the 
electromagnetic vector potential further increase the difference between the action (or accumulated phase) of electrons 
and holes, it will increase the degree of imperfectness in the retro-reflection. 

The present work, and in particular approximation Scheme B, provides insight into the general question of when 
electron dynamics should be handled separately from the hole dynamics in inhomogeneous superconductors. In doing 
this, it also provides a semiclassical framework for studying the effects of electron/hole symmetry breaking beyond 
Andreev approximation. From the point of view of quantum chaos, such electron/hole differences lead to the novel 
dynamics featuring in the present work. Phenomena associated with this should be accessible via experiment and, 
indeed, the extremely recent Regenbsurg experiments are paving the way to a thorough experimental exploration of 
the novel quantal dynamics of electrons and holes in Andreev billiards. 
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APPENDIX A: APPLICATION OF BOUNDARY INTEGRAL TECHNIQUE TO THE HELMHOLTZ AND 

BDG WAVE EQUATIONS 

The underlying strategy employed in this Paper is the boundary integral technique EJ, the origin of which was 
Fredholm's analysis of the existence of soluions of the interior Laplace problem subject to Dirichlet boundary condi- 
tions. In this scheme, Fredholm transformed the task of solving the Laplace partial differential equation (subject to 
Dirichlet boundary conditions) to one of solving a certain integral equation residing on the boundary. This prompted 
Fredholm to develop the theory of what are now known as Fredholm integral equations and, in particular, to prove 
the existence of a solution of the corresponding Laplace problem. 

In the present context of spectral geometry, the virtue of this boundary integral technique is that it allows one 
to harness the piecewise homogeneity of the system (and the corresponding simplicity of the fundamental Green 
functions in the locally homogeneous regions) and, thereby, to study the physical implications of the boundary in as 
direct and a natural manner as possible. 

The aim of this appendix is to provide a guide to the boundary integral technique, beginning with the simplest 
setting and working towards the setting of the BDG eigenproblem. When discussing the simplest settings we shall 
borrow heavily from Refs. |^2|,Q. The elaborations that we shall be needing for the BDG setting arise from (i) the 
multicomponent nature of the eigenproblem, and (ii) the presence of matching rather than boundary conditions. We 
mention that Jackson Q| gives a highly readable discussion of the physics of the potential discontinuities that are a 
pivotal feature of boundary integral techniques. 

1. Review of elementary ingredients 

In this section we shall discuss the origin of the discontinuities in the three-dimensional potential and the fields 
generated by surface charge (which we call single-layer) and dipole (which we call double-layer) densities, as well as 
present derivations of explicit formulas quantifying such discontinuities. Such formulas will become useful in the next 
section, where we discuss parametrizations of wave functions in terms of these single and double layers. 
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Before turning to the derivation, we consider a simple example which contains the essential features: a planar 
charge layer with constant charge density v. Without loss of generality, let us assume that the charge layer lies in the 
xy plane, so that the normal direction nis z. Then the potential <p and the field E are given by the following surface 
integrals: 

47T<,5(x) = v [ da a -. r, (Ala) 

J |x-a| 

4 7 rE(x) = -,|^V,( R ^) =vjda a ^^, (Alb) 

where x = (x, y, z), a = (x 1 ' , y', 0), and da a = dx' dy' . Owing to translational invariance in the xy plane, y(x) = ip(z) 
and E(x) = zE(z). Without loss of generality, let us choose x = y = 0, and focus on E(z). By introducing the polar 
coordinates {x',y') = (r cos 6*,r sin #), the integral for E(z) becomes 

/ Z/ f°° Z V 

rdrdd Mz* +r>)S/z = 2j rdr (z2 + Z r2 )3/2 - r gn{z) (A2) 

For z > 0, E(z) is independent of z and equal to v/2; for z < 0, E{z) is independent of z, however it is equal to 
—v/2 and for z = 0, E(z) vanishes. The mathematical origin of this discontinuity is the noncommutativity of the 
limit z — > and the surface integral in Eq. (A2). Thus the value of ip is continuous as it approaches the surface but 
its normal derivative (in this case partial derivative with respect to z) is not. This discontinuity can be summarized 
by the equation 

lim cp = (p\g= , lim V 5 = (p\g=o, (A3) 

2^0+ z^Q- 



lim ftt = -- + ^L 

z^o+ dz 2 dz 



dip v dp 
z^o- oz 2 oz 



(A4) 



Now consider a generic charge layer v(a) on the surface dV, and focus on the potential generated by via the 
Helmholtz (instead of Coulomb) Green function 

<p(x)= / dcr Q 5 H (x, a) v(a). (A5) 

JdV 

As in the case of the simple example of homogeneous planar charge layer, the potential generated by this generic 
charge layer is continuous, 



_lim <p(x) = / da a g H (/3, a) v{a) 



(A6) 



but its normal derivative is not. In order to see this, consider the normal derivative • </?(x) as x in V tends to a 
generic point on dV, which we denote by (3, along the interior normal and divide the domain of the surface integration 
into two parts: (i) a small region D$ = C$ f] dV, where Cs is a sphere of radius S around (3 and (ii) remaining domain 
V s = dV -T>s, then 

[ d f d 

_lim nn ■ Vj, pipe) = lim / da a — g n {f3 + e ng, a) v(a) + lim / da a — g U {(3 + e n«, ex) v(a) (A7) 

where e is the perpendicular distance between x and (3. The virtue of this separ atio n is that the singularity of g H (/3, a) 
at f3 = a is now contained in the first term on the right hand side of the Eq. (A7). As the integrand of second term 
on the right hand side of this equation is free of singularities, the limit can be taken inside the integral sign. Moreover 
within T>$ and for S and e very small, 

^ff H (/3 + e n ff , Q )„ 4 ^ 2 ^ e2)3/2 , p^|«-/3| 2 . (A8) 

Then the integral over the domain T>s can be evaluated for small 5 as follows 

f d f 2lT f s el 

lim dcr a —g ii {P + enp,a)v{a) = limi/(/3) ^ dd ^ P d P ^ p2 + £ 2)3/2 = (A9) 
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Now that we have evaluated the singular part of the surface integral we take the limit 5 tends to zero: 

Jim n • V* tp(x) = -1/09) + lim / da a • V g H {/3, a) v{a) (A10) 
= \v(P)+ I da a n fj - V p g H ((3, a) iy(a). (All) 

1 JdV 

We note that this discontinuity originates again from the noncommutativity of limit x — + /3 and the surface integral. 
Moreover the calculation of the amount of this discontinuity involves only the form of <? H for small values of distance 
between its arguments and this form is simply the Coulomb Green function. Thus, in essence, this discontinuity is 
the same as the discontinuity in the simple case of a planar, electrostatic-charge layer discussed at the beginning of 
this section. 

Let us now consider the potential generated by a generic dipole density n(ot) on dV, viz. 

ip(x) = / da a n a ■ V Q g H (x, a) fj,(a). (A12) 

JdV 

Notice the similarity of this form to that of the normal derivative of the potential generated by the charge layer, only 
difference being the normal derivative acting on the second rather than the first argument. However, as the G H is 
a function of the difference between its arguments, one expects a similar discontinuity on the potential across the 
surface. Let us now consider the case in which x approaches to a surface point denoted by /3 from V. Indeed, 

_lim tp(x) = _lim / d(7 a n 'V a (/ H (x, a)/j(a) (A13) 
v^pedv xev^peov Jdv 

_lim / da a n Q • V a g H (x, a) fi(a) + _lim / da a n a • V a <? H (x, a) /x(a) (A14) 
v^pedvJvs xev^pedvJv s 

lim / da a • V x <? H (x, a) /x(a) + _lim / da a n a ■ V Q 5 H (x, a) fi(a), (A15) 
^pedvJ-Ds xev-tpedvJvs 



xt 



where in order to get to the third line we have used V x f/ H (x, x') = — V^g^x, x'), and that in T>s, as 5 goes to zero, 
n a — * rip. Notice that the first term in the right hand side of Eq. ( A15) is equal to the first term in Eq. (A9). Thus 
in the limit S goes to zero, we have: 

Jim (p(x) = --/i(/3) + / da a n a ■ V a g H ((3 , a) n(a) . (A16) 
xev^p&ov & Jdv 

The discontinuity in tp as x approaches to /3 from V (rather than V), is obtained similarly, the result is: 

lim p(x) = -fi(P) + / da a n a ■ V a g H (f3 , a) fi(a) . (A17) 
xev^pedv z J dv 



2. Single- and double-layer parametrizations of wave functions; Jump conditions 

The first ingredient needed for the construction of a MSE is the parametrization of a wave function in terms of 
single or double layers. As this aspect of classical potential theory might be unfamiliar to some readers, we first 
illustrate it in the simpler setting of the Helmholtz wave equation, before turning to the BDG wave equation. 



a. Parametrizations for one- component Helmholtz wave functions 
Consider the Helmholtz wave equation, 

(V 2 + E) <p{x) = 0, (A18) 

for the wave function y(x) in the region V (which we define to be the region outside some region V bounded by the 
closed surface dV). Then, from potential theory jf2|,fl3| it is known that the solutions ip{x) can be parametrized in 
terms of a function v{a) defined only on dV, via the integral 
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<y9( X ) = 



/ do- a .g H (x- a)v(a). (A19) 

JdV 



Here, x and x' are positions lying in V, Greek letters, such as a, represent vectors on the boundary dV (as they do 
throughout this Paper), and g H (x — x') is the fundamental Green function for the Helmholtz wave equation, which 
satisfies 

(V 2 +£) g H (x-x') =5( 3 )(x-x'). (A20) 

One can interpret the parametrization by saying that the wave function y(x) is the Helmholtz wave function due to 
a single layer of charge of surface density v{ot). 

Where does this parametrization come from? First, note that the wave function ( Al9| ) does indeed satisfy Eq. ( Al§| ). 



To s ee th is, observe that x and a are never coinc ident (a lying on dV but x lying in V) so that <? H (x — a) solves 
Eq. (|Al|) for any a, and the parametrization |A19| is simply a superposition of such solutions. Second, recall that by 



Green's theorem one has 

ip(x) 



/ da a d a g A (x,a)ip(a) - da a g A (x,a) d a ip(a), (A21) 

JdV JdV 



where g A (x, a) is any Hclmhotz Green function (i.e. not necessarily the fundamental one). In the most common 
setting, one then chooses the Green function that satisfies the homogeneous version of the boundary condition on 
ip, thus eliminating all absent boundary information and arriving at an expression for </?(x) in terms of information 
known about ip on dV. Here, instead, one takes a different tack. One selects for g A the fundamental Green function 



g , jettisons the first contribution to the right hand side of Eq. (A21), and accommodates for this by replacing the 
boundary information on dip by the (as-yet unknown) single-layer v(a). 

Representations of wave functions by surface integrals are available in other settings, too. We have considered wave 
functions satisfying the Helmholtz wave equation outside the region V (i.e. the so-called exterior problem). One can 
also consider wave functions satisfying the Helmholtz wave equation inside the region V (i.e. the so-called interior 
problem). 

Furthermore, one can parametrize wave functions satisfying the Helmholtz wave equation in other ways. For 
example, consider wave functions inside the region V, which one can parametrize as 



tp(x) 



/ do- a d a g H (x- a)n(a). (A22) 

JdV 



In this case ip(x) is the Helmholtz wave function due to a layer of dipoles on dV of local strength (J,(et) and local 
orientation normal to dV at each point a. (We denote such surface normal vectors as n Q , and adopt the convention 
that they point towards the interior of V.) One can, of course, regard this dipole layer as consisting of two single 
layers, vanishingly close to one another and locally carrying opposite charges, in the limit that the charges become 
large and the layer separation becomes correspondingly small. Such layers are referred to as a double layers. The 
normal component of the gradient acting on the Green function accounts for the fact that this parametrization fea ture s 
opposing, van ishin gly close, layers. As with the single-layer parametrization (A19), that the para metriz ation ( A22j ) 



satisfies Eq. ( A18 ) follows because x an d a are never coincident, so that g (x — a) solves Eq. ( Ala ) for any a 



Motivation for the parametrization (A22) also follows from consideration of Green's theorem, Eq. (A21), but with the 
second contribution on the right hand side being jettisoned, rather than the first. 

There, of course, remains the issue that whether all solutions can be expressed in terms of these parametriza- 
tions. It turns out that if one uses the double-layer parametrization for interior wave functions, and the single-layer 
parametrization for exterior wave functions then any solution can be thus parametrized. The converse problem 
(i.e. parametrizing interior wave functions using double layers and exterior wave functions using single layers) is also 
possible, provided certain supplementary conditions are satisfied; see, e.g., Refs. fl^Jfll 

We note that the strategy that we are adopting can be implemented in more general settings. For example, one 
might consider the case of disconnected superconducting regions connected by normal regions, and thus address the 
issue of Josephson tunneling between them. One might also consider disconnected normal regions connected by 
superconducting regions, and thus address the issue of single-particle tunneling between them. 

The utility of these parametrizations is that they can be used to transform partial differential equations in V or 
V into Fredholm integral equations that reside on dV, and as we shall see in the following section, such integral 
equations prove useful in some cases, especially if the integral equation is of the second type, i.e., if an iterative 
solution is possible. 

We now turn to the second ingredient needed for the construction of a MSE, viz., jump conditions. It is an important 
result of potential theory that the parametrization of wave functions in terms of single and double layers, such as those 
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, leads to representations of wave functions that behave in a singular fashion for field points 
x on the surface dV. It is precisely this singular behavior, and the attendant jump conditions, that are responsible 
for the utility of these parametrizations and, as we shall see shortly, lead to the formulation of integral equations for 
the layer strengths, single or double, known as boundary integral equations. These equations incorporate what ever 
boundary conditions one wishes to impose on the wave function. By solving boundary integral equations one arrives 
at layer strengths that parametrize the wave functions. 

There are three virtues to this boundary integral equations formulation. First, the boundary integral equations 
reside solely on the boundary dV (i.e. on a manifold of dimension one fewer than the original wave equation). From 
the computational standpoint it is economical to formulate a problem in terms of functions that reside on lower- 
dimensional manifolds (and hence depend on fewer variables). Second, from the theoretical standpoint, existence 
theorems have been established for broad classes of integral equations, often encompassing the boundary integral 
equations that emerge from specific examples. (In fact, it was the goal of establishing exitence theorems for solutions 
to the Laplace equation in various settings — interior or exterior, Dirichlet or Neumann boundary conditions — that 
inspired Frcdholm to develop the boundary intergal equation approach to potential theory, and subsequently to 
develop the theory of what we now know as Fredholm integral equations.) And third, from the physical standpoint, 
boundary integral equations and their iterative solution allow one to organize the computation of wave functions in 
terms of the multiple scattering of waves from interfaces that separate spatially homogeneous regions, along with free 
propagation between those scattering events. Thus, one is in a position to focus on the boundary scattering events, 
and thereby to focus on the geometry of the boundary and the implications of its shape for the physical problem at 
hand. In essence, we are invoking the piecewise homogeneity of the system to "integrate up" our description of it, 
leaving us with the need to consider one fewer independent variable. As a result of this "integrating up" we depart 
from a purely local description, in terms of differential equations, and arrive at a nonlocal formulation in terms of 
integral equations. 



given in Eqs. (|AiqjA22| 



b. Jump conditions for one-component Helmholtz wave functions 



Consider the single-layer parametrization for Helmholtz wave functions in the region V, given by Eq. (A19). Now, 
it is known from potential theory that this parametrization is singular as x goes to any value (3 on the boundary (see 
App. AS). Specifically, ip(x) is continuous whereas • V x ip(x) is not, but the discontinuity of the latter has a known 
and useful form: 



xev- 



lim / da a <? H (x — a) v{ol) = I da a g H (/3 — a) v(a); (A23a) 
•->0£dvj gv J dv 

lim np-W x da a g H (x-a)v(a.) = lv{l3)+ da a n -V >S H (/3 - a)v{a). (A23b) 
■^pedv J gv I J gv 



Now consider the double-layer parametrization of the Helmholtz wave functions on the region V, given by Eq. (|A22j ). 
It is known from p oten tial theory that this parametrization is also singular as x goes to any value of (3 on the 
boundary (see App. |A 2| ). However, in this case </?( x ) itself is discontinuous, and rip ■ V x </?(x) is even more singular, 
the discontinuity of ip{x) being given by 



xev- 



lim / da a d a g H (x-a)n(a) = l-n{(3)+ da a g H ((3 - a) fi(a). (A24) 



Although rip ■ V x (p(yi) diverges on dV, it is yet a further result from potential theory that the values of the limits 
of this quantitity, as x approaches any point (3 on dV either from V or from V, are equal to one another. It will, 
therefore, prove to be convenient to redefine the quantity rip ■ V x <y9(x)| x=/ g as its limiting value: 

/ da a dt d a g ii (f3 — a)u(a) = lim rip-V x / da a d a g H (/3 — a) v{ot). (A25) 

JdV x£V or V^/3£dV J dV 

Via this definition the normal derivative • 'V x tp(x) is continuous across dV. 



c. Parametrizations for two-component BDG wave functions 

We now turn to the issue of parametrizing solutions of the BDG equation (^) in terms of single and double layers, 
bearing in mind the two-component nature of the wave functions. For the sake of concreteness, as well as experimental 
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relevance, we focus on the setting of an Andreev billiard, so that space is partitioned into two regions by the surface 
dV, with A(x) = in the region V inside dV (i.e. the normal-metal region) and A(x) = Ao in the region V outside 
dV (i.e., the superconducting region). However, the are no obstacles of principle in applying the present techniques 
in other settings, provided they comprise regions of space in which A(x) is constant, separated by surfaces on which 
A(x) varies discontinuously [i.e. A(x) must be piecewise constant]. The present techniques may also be applied in 
settings in which other physical parameters vary in a piecewise continuous fashion. 

First, consider the normal-state interior of an Andreev billiard [i.e. a region V surrounded by a surface dV in which 
A(x) = 0] in which the BDG wave functions satisfy 

n \ /iiYvlN /n\ 

(A26) 

Solutions of this equation can be parametrized in terms of the two-component double layer [J,(et) in the following 
way: 

* N (x)=/ da a d a G N (x-a) -n(a). (A27a) 

JdV 

Here, G N (x — x') is the fundamental Green function for the BDG wave equation in the absence of a pair potential, 




and is given explicitly by Eq. (2.19). Such solutions could also be parametrized in terms of the two-component single 
layer u{a), as 

* N (x) = / da a G N (x - a) ■ v(a), (A27b) 

JdV 



but it will prove convenient to adopt the former parametrization, Eq. (A27a), rather than the latter, Eq. ( |A27b ). 



The two-component layers, (a(o.) and u (a), reflect the two-component (i.e. electron and hole) nature of the wave 
functions. If one were solely concerned with the case of a normal region, this two-component description would be 
redundant: as no matrix elements of the Hamiltonian would connect upper and lower components there would be no 
need to adopt a language that embraces wave functions that describe coherent superpositions of electron and hole 
states. However, in an Andreev billiard the normal region is surrounded by a superconductor, the pair potential of 
which provides precisely the matrix element connecting electron and hole wave functions. Therefore it is necessary to 
adopt this two-component language for the normal region. 

In the superconducting exterior of the Andreev billiard (i.e. the region V outside the surface dV in which A(x) = Ao) 
the BDG wave functions satisfy 



' V 2 - fi + E -iA 

iA a -V 2 + /i + E 




(A28) 



Solutions of this equation can be parametrized in terms of of the two-component single layer v(a) in the following 
way: 

* s (x)^ = J^da a G s (x- <*)■„(<*), (A29a) 

where G s (x — x') is the fundamental Green funct ion for the BDG wave equation in the presence of a homogeneous 



pair potential Ao, and is given explicitly by Eq. (2.27b). Two-component layers are mandatory here, inasmuch as 
each component of the wave function is determined by both components of a layer, owing to the presence of the 
pair-potential. Such solutions could also be parametrized in terms of the two-component double layer /x(ck), as 

* s (x) = H*J J da a 8 a G s (x - a) ■ M (a), (A29b) 

except that under certain circumstances the parametrization must be augmented by an additional term; see, e.g., 
Ref. E3[. However, it will prove adequate for us to stick with the former parametrization. 
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d. Jump conditions for two-component BDG wave functions 



We now echo for the case of BDG wave functions the discussion, given in App. A 2 t , of the behavior of single- and 
double-layer parametrizations of Helmholtz wave functions in the vicinity of the su rface dV. Consider the double-layer 
parametrization of the BDG wave functions on the region V given by Eq. ( A27a ). This parametrization is singular, 
as x goes to any point (3 on the boundary (see App. A 2), inasmuch as <fr N (x) is discontinuous, and rip • V 2; ^ ,N (x) 
is even more singular. From the form of G N (x — x'), Eq. ( 2.21 ), it can be shown that the discontinuity of $ N (x) is 
given by 



lim f d<r a d a G N (x — a) ■ fJ.(a) = —cr^ ■ (i(f3) + [ da a d a G N (/3 — a) ■ ju(ct) 
xev^pedv J dv 2 J dv 



(A30a) 



Although xip ■ V x $ w (x) diverges on dV, it can also be shown that that the values of the limits of this quantitity, as 
x approaches any point (3 on dV either from V or from V, are equal to one another. It will, therefore, prove to be 
convenient to redefine the quantity ■ V x $ N (x)| x= ^g as its limiting value: 



/ da a d t d a G N (/3 — a) ■ n(a) = lim rig • V x / da a d a G N (x — a) ■ fx(a). 

JdV xeV or v^pedv JdV 



(A30b) 



Now consider the single-layer parametrization of the BDG wave functions in the region V given by Eq. ( A27b). This 
parametrization is singular, as x 6 V goes to any value of (3 € dV (see App. A 2). Specifically, 3? S (x) is continuous 



wher eas rig • (x) is not, but the discontinuity has a known and useful form. Indeed, from the form of G (x — x'), 



Eq. ( 2.27b| ), it can be shown that 



_lim / da a G s (x — a) ■ v(a) 

xgv^/3g9v7ov 



da a G s (f3 - a) ■ v{a); 



t)V 



Jim np-V x / da a G (x - a) ■ u(a) = -<r 3 • v{0) 
xGV^/3eav Jdv 1 



da a np ■ Vp G s (/3 — a) • u(a). 



(A31a) 
(A31b) 



The important point here is that all discontinuities of these parametrizations are generated solely by components 
of G N,S proportional to cr^. Other components of G N ' S are proportional to the scalar Green function composition 

N S N S 

(<7_l_' — gJ ), and whatever discontinuity might be generated by the + term is cancelled by a corresponding one 
gener ated b y the — term . Therefore the discontinuities generated by the parametrizations have the forms given in 
Eqs. (|A30a|) and (lAfTg). 



APPENDIX B: (d - l)-DIMENSIONAL FOURIER TRANSFORMS 

In this appendix we introduce the d — 1 dimensional Fourier transform for functions whose arguments are d di- 
mensional vectors, which we use extensively throughout the text. Mainly, our functions of interest will be the Green 
functions g± and the functions related to them, such as dg± . 

The d—1 dimensional Fourier transform (and its inverse) of a function /(x) of a d dimensional vector x, are defined 

by 



f(K,z)= / dx // e- x ///(x), (Bla) 
Jv 

/(x) = (2vr)- d+1 J dne- iK ^ f(n,z), (Bib) 

where (x //, z) = x, and x // and k are vectors on the d—1 dimensional hyperplane V perpendicular to the z axis. We 
now evaluate the (d — l)-dimensional Fourier transforms of the functions that are used throughout this Paper. We 
begin with the Helmholtz Green function g(x;k 2 ): 
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where a(n) = V 'k 2 — k 2 such that Re(a(re)) > 0. The evaluation of the (d — l)-dimensional Fourier transform of 
dSg(x; k 2 ) is very similar, except that now x = (x//,0) (i.e. both arguments of the Green function lie on V and the 
normal direction is z): 



d5g(n; k 2 ) 



dx //e ^d5g((x //7 0);k 2 ) 



dx// e lK ' x // 



dz 2 



dpz pi e 



2 p ip z z 



dp e ip //' x // +iPiZ 
{2n) d p 2 - k 2 



2ir p 2 z + a(n) 2 



z=0 



(B3) 



APPENDIX C: PROPAGATION OUTSIDE THE N REGION AND NONCONVEX SHAPES: 

CANCELLATIONS 



In the present section we study cancellations between terms in the MSE, which occur, in the large kpL limit, when 
one (or more) of the homogeneous region(s) constituting the billiard are nonconvex. In the sample billiard shown in 
Fig. m the S region is nonconvex, whereas the N region is convex. Another example in which such cancellations arise 
is provided by antidot-billiard geometries, in which S regions are embedded in an N region. For such cases, the N 
region is certainly nonconvex, and so may be the S region. In the former case (i.e. convex N; nonconve x S), t hese 
cancellations eliminate terms that include finite-outside-propagation, validating the claim made in Sec. IV D that 
periodic orbits that include such paths do not contribute at leading order to the semiclassical DOS oscillations. 



a 




FIG. 15. Example of a nonconvex inside and outside regions. Also shown is a propagation directly from a to (3 and 
two-reflection correction. 



In order to bring to the fore the physics underlying these cancellations, consider a much simpler case: the MRE for 
the Helmholtz Green function satisfying homogeneous Dirichlet boundary conditions Jf5[ . By applying the methods 
described here, one can express this Green function in terms of the following MRE: 



. 9 D (x,x')=ff H (x,x')-2 f % H (x,a) 5 H (a,x')+4 f 9.g H (x, 0) % H (/3, a) g H (a, x' 

JdV JdV 



(CI) 



where <? H is the homogeneous Helmholtz Green function. If the domain V over which g° is defined is nonconvex then 
one may be concerned by the presence of amplitudes involving propagations between pairs of boundary points for 
which all or part of this propagation lies outside of V (see Fig. Hq). On physical grounds, one expects that in the large 
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kL limit the obstacle between these points would suppress such amplitudes. Balian and Bloch showed that this is 
indeed the case. To see this, we follow Balian and Bloch and consider g° for a planar boundary. Then the MRE (CI) 
terminates at the second term: 



. 9 D (x,x')=g H (x,x')-2 f d g u (x >a )g u ( a , X '). 

JdV 



(C2) 



On the other hand, by using the method of images one has 

.g D (x,x')=.9 H (x,x')-.9 H (x,x'), 

where x' is the mirror image of x' with respect to the planar boundary. Comparison of Eq. ( |C3| ) with Eq. ( |C2| ) 
produces the relation 



(C3) 



5 H (x,x')=2/ 9. 9 H (x,a) 5 H (a,x'). 

JdV 



(CM) 



As the length \a — x'| equals the length \a — x'|, one can replace x' by x' in relation (C4) and obtain a second relation 



H (x,x')=2/ 9.9 H (x,a) 5 H (a,x'). 

JdV 



(C5) 



Now suppose dV is not planar. Although relations ( |C4| ) and (C5) no longer hold, in the large kL limit the dominant 
contribution comes from the stationary-phase point and, thus, Eq. (|C5| ) becomes exact as kL tends to infinity. That 
Eq. (C4) is not exact in this limit is due to the fact that the fluctuation determinant in this case (viz. when x and x' 
are on the same side of the surface) depends on the curvature of the surface at the stationary-phase point. 

Now let us return to what happens for nonconvex surfaces. In particular, consider a term in MRE in which part of 
the amplitude has propagation between boundary points a. and f3, where part of the line joining a and (3 lies outside 
V and intersects dV at the points 70 and <5o (see Fig. |l5[). In relation to this term, consider three related terms with 
higher numbers of reflections: (i) the term with one more reflection near 70, (ii) the term with one more reflection 
near Sq, and (iii) the term with two more reflections near 70 and Sq. In the MRE, all four are in the sum. The sum 
of the part of the amplitude involving direct propagation from a to (3 and the one with one more reflection near 70 is 



dg H ( a ,(3)-2[ 3 5 H (a,7)d 5 H (7,/3). 

JdV 



(C6) 



By virtue of Eq. (C5) this sum vanishes in the kL — ► 00 limit. The same holds for the sum of remaining two terms. 
Thus in the limit kL — > 00 sum of all terms involving paths that lie partially outside of V vanishes. 

These considerations extend readily to Andreev billiards. In order to see this, we must extend the identity (|C5| ). 
This can be achieved as follows: 



2 / dG (x, a) <r 3 G (a, x) = 2 

lav Jdv 



<9ff+(x,a).9 + (a,x) 

-dg-(x,a)g-(oc,5i) 



_ ( ff+(x,x) 






-5-(x,x) 



= G N (x,x), 



(C7) 



where, as in the case of Hclmholtz Green functions, dV is a planar surface and x and x lie on different sides of dV, 
and in going to the second line we have made use of Eq. (C5). A similar identity holds for G s 



G b (x,x)=2/ <9G b (x,a)er 3 G b (c*,x) = 2 / G b (x, a) cr 3 SG b {a, x 



av 



(C8) 



av 



This identity can be obtained by performing the matrix mu ltiplic atio n and using Eq. ( |C5| ) in each matrix element. 
As in the case of Helmholtz Green functions, identities (C7) and ( |CS| ) hold in the large k-pL limit, even if dV is not 
planar. As we did for the Helmholtz case, consider a term in MSE in which part of the amplitude has propagation 
between boundary points a. and (3, where part of the line joining a and (3 lies outside V and intersects dV at the 
points 70 and Sq. However, now we have additional terms owing to the possibilit y of propagation between any two 



boundary points involving either G or G . By using the identities (C7) and (C£) it is not hard to see that, in the 



large kpL limit, a propagation from a to (3 is not cancelled only if the line joining these points lies totally in either V 
or V, with the propagation involving the homogeneous Green function appropriate to the corresponding region. Thus, 
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in the semiclassical limit the terms that survive in the MSE consist of pure reflection or pure transmission/tunneling. 
For the example at hand, the surviving term would be the one with normal propagation from a to 7, superconducting 
propagation from 7 to S (i.e. tunneling), and then normal propagation from d to (3. 
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